{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kelvin-Helmholtz Instability\n",
    "\n",
    "<img src=\"http://www.flame.org/~cdoswell/SuptorRoles/Vortex%20Sheet%20Model%20copy.jpg\" width=\"254\" height=\"300\" />\n",
    "\n",
    "(image: Chuck Doswell)\n",
    "\n",
    "We will simulate the incompressible Kelvin-Helmholtz problem.  We non-dimensionalize the problem by taking the box height to be one and the jump in velocity to be one.  Then the Reynolds number is given by\n",
    "\n",
    "$$ \\mathrm{Re} = \\frac{U H}{\\nu} = \\frac{1}{\\nu}. $$\n",
    "\n",
    "We use no slip boundary conditions, and a box with aspect ratio $L/H=2$.  The initial velocity profile is given by a hyperbolic tangent, and only a single mode is initially excited.  We will also track a passive scalar which will help us visualize the instability.\n",
    "\n",
    "First, we import the necessary modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import h5py\n",
    "import time\n",
    "from IPython import display\n",
    "\n",
    "from dedalus import public as de\n",
    "from dedalus.extras import flow_tools\n",
    "\n",
    "import logging   \n",
    "logger = logging.getLogger(__name__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To perform an initial value problem (IVP) in Dedalus, you need three things:\n",
    "\n",
    "1. A domain to solve the problem on\n",
    "2. Equations to solve\n",
    "3. A timestepping scheme\n",
    "\n",
    "## Problem Domain\n",
    "\n",
    "First, we will specify the domain.  Domains are built by taking the direct product of bases.  Here we are running a 2D simulation, so we will define $x$ and $y$ bases.  From these, we build the domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Aspect ratio 2\n",
    "Lx, Ly = (2., 1.)\n",
    "nx, ny = (192, 96)\n",
    "\n",
    "# Create bases and domain\n",
    "x_basis = de.Fourier('x', nx, interval=(0, Lx), dealias=3/2)\n",
    "y_basis = de.Chebyshev('y',ny, interval=(-Ly/2, Ly/2), dealias=3/2)\n",
    "domain = de.Domain([x_basis, y_basis], grid_dtype=np.float64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last basis ($y$ direction) is represented in Chebyshev polynomials.  This will allow us to apply interesting boundary conditions in the $y$ direction.  We call the other directions (in this case just $x$) the \"horizontal\" directions.  The horizontal directions must be \"easy\" in the sense that taking derivatives cannot couple different horizontal modes.  Right now, we have Fourier and Sin/Cos series implemented for the horizontal directions, and are working on implementing spherical harmonics.\n",
    "\n",
    "## Equations\n",
    "\n",
    "Next we will define the equations that will be solved on this domain.  The equations are\n",
    "\n",
    "$$ \\partial_t u + \\boldsymbol{u}\\boldsymbol{\\cdot}\\boldsymbol{\\nabla} u + \\partial_x p = \\frac{1}{{\\rm Re}} \\nabla^2 u $$\n",
    "$$ \\partial_t v + \\boldsymbol{u}\\boldsymbol{\\cdot}\\boldsymbol{\\nabla} v + \\partial_y p = \\frac{1}{{\\rm Re}} \\nabla^2 v $$\n",
    "$$ \\boldsymbol{\\nabla}\\boldsymbol{\\cdot}\\boldsymbol{u} = 0 $$\n",
    "$$ \\partial_t s + \\boldsymbol{u}\\boldsymbol{\\cdot}\\boldsymbol{\\nabla} s = \\frac{1}{{\\rm Re}{\\rm Sc}} \\nabla^2 s $$\n",
    "\n",
    "The equations are written such that the left-hand side (LHS) is treated implicitly, and the right-hand side (RHS) is treated explicitly.  The LHS is limited to only linear terms, though linear terms can also be placed on the RHS.  Since $y$ is our special direction in this example, we also restrict the LHS to be at most first order in derivatives with respect to $y$.\n",
    "\n",
    "We also set the parameters, the Reynolds and Schmidt numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "Reynolds = 1e4\n",
    "Schmidt = 1.\n",
    "\n",
    "problem = de.IVP(domain, variables=['p','u','v','uy','vy','S','Sy'])\n",
    "problem.parameters['Re'] = Reynolds\n",
    "problem.parameters['Sc'] = Schmidt\n",
    "problem.add_equation(\"dt(u) + dx(p) - 1/Re*(dx(dx(u)) + dy(uy)) = - u*dx(u) - v*uy\")\n",
    "problem.add_equation(\"dt(v) + dy(p) - 1/Re*(dx(dx(v)) + dy(vy)) = - u*dx(v) - v*vy\")\n",
    "problem.add_equation(\"dx(u) + vy = 0\")\n",
    "problem.add_equation(\"dt(S) - 1/(Re*Sc)*(dx(dx(S)) + dy(Sy)) = - u*dx(S) - v*Sy\")\n",
    "problem.add_equation(\"uy - dy(u) = 0\")\n",
    "problem.add_equation(\"vy - dy(v) = 0\")\n",
    "problem.add_equation(\"Sy - dy(S) = 0\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we are using this first-order formalism, we define auxiliary variables `uy`, `vy`, and `Sy` to be the $y$-derivative of `u`, `v`, and `S` respectively.\n",
    "\n",
    "Next, we set our boundary conditions.  \"Left\" boundary conditions are applied at $y=-Ly/2$ and \"right\" boundary conditions are applied at $y=+Ly/2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "problem.add_bc(\"left(u) = 0.5\")\n",
    "problem.add_bc(\"right(u) = -0.5\")\n",
    "problem.add_bc(\"left(v) = 0\")\n",
    "problem.add_bc(\"right(v) = 0\", condition=\"(nx != 0)\")\n",
    "problem.add_bc(\"left(p) = 0\", condition=\"(nx == 0)\")\n",
    "problem.add_bc(\"left(S) = 0\")\n",
    "problem.add_bc(\"right(S) = 1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we have a special boundary condition for the $k_x=0$ mode (singled out by `condition=\"(nx==0)\"`).  This is because the continuity equation implies $\\partial_y v=0$ if $k_x=0$; thus, $v=0$ on the top and bottom are redundant boundary conditions.  We replace one of these with a gauge choice for the pressure.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Timestepping\n",
    "\n",
    "We have implemented a variety of multi-step and Runge-Kutta implicit-explicit timesteppers in Dedalus.  The available options can be seen in the [timesteppers.py module](https://github.com/DedalusProject/dedalus/blob/master/dedalus/core/timesteppers.py).  For this problem, we will use a third-order, four-stage Runge-Kutta integrator.  Changing the timestepping algorithm is as easy as changing one line of code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = de.timesteppers.RK443"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial Value Problem\n",
    "\n",
    "We now have the three ingredients necessary to set up our IVP:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2020-10-30 21:55:41,722 pencil 0/1 INFO :: Building pencil matrix 1/96 (~1%) Elapsed: 0s, Remaining: 5s, Rate: 1.8e+01/s\n",
      "2020-10-30 21:55:42,027 pencil 0/1 INFO :: Building pencil matrix 10/96 (~10%) Elapsed: 0s, Remaining: 3s, Rate: 2.8e+01/s\n",
      "2020-10-30 21:55:42,271 pencil 0/1 INFO :: Building pencil matrix 20/96 (~21%) Elapsed: 1s, Remaining: 2s, Rate: 3.3e+01/s\n",
      "2020-10-30 21:55:42,489 pencil 0/1 INFO :: Building pencil matrix 30/96 (~31%) Elapsed: 1s, Remaining: 2s, Rate: 3.6e+01/s\n",
      "2020-10-30 21:55:42,654 pencil 0/1 INFO :: Building pencil matrix 40/96 (~42%) Elapsed: 1s, Remaining: 1s, Rate: 4.0e+01/s\n",
      "2020-10-30 21:55:42,842 pencil 0/1 INFO :: Building pencil matrix 50/96 (~52%) Elapsed: 1s, Remaining: 1s, Rate: 4.3e+01/s\n",
      "2020-10-30 21:55:43,020 pencil 0/1 INFO :: Building pencil matrix 60/96 (~62%) Elapsed: 1s, Remaining: 1s, Rate: 4.4e+01/s\n",
      "2020-10-30 21:55:43,187 pencil 0/1 INFO :: Building pencil matrix 70/96 (~73%) Elapsed: 2s, Remaining: 1s, Rate: 4.6e+01/s\n",
      "2020-10-30 21:55:43,397 pencil 0/1 INFO :: Building pencil matrix 80/96 (~83%) Elapsed: 2s, Remaining: 0s, Rate: 4.6e+01/s\n",
      "2020-10-30 21:55:43,565 pencil 0/1 INFO :: Building pencil matrix 90/96 (~94%) Elapsed: 2s, Remaining: 0s, Rate: 4.7e+01/s\n",
      "2020-10-30 21:55:43,670 pencil 0/1 INFO :: Building pencil matrix 96/96 (~100%) Elapsed: 2s, Remaining: 0s, Rate: 4.8e+01/s\n"
     ]
    }
   ],
   "source": [
    "solver =  problem.build_solver(ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we set our initial conditions.  We set the horizontal velocity and scalar field to tanh profiles, and using a single-mode initial perturbation in $v$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Field 140447541869648>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = domain.grid(0)\n",
    "y = domain.grid(1)\n",
    "u = solver.state['u']\n",
    "uy = solver.state['uy']\n",
    "v = solver.state['v']\n",
    "vy = solver.state['vy']\n",
    "S = solver.state['S']\n",
    "Sy = solver.state['Sy']\n",
    "\n",
    "a = 0.05\n",
    "sigma = 0.2\n",
    "flow = -0.5\n",
    "amp = -0.2\n",
    "u['g'] = flow*np.tanh(y/a)\n",
    "v['g'] = amp*np.sin(2.0*np.pi*x/Lx)*np.exp(-(y*y)/(sigma*sigma))\n",
    "S['g'] = 0.5*(1+np.tanh(y/a))\n",
    "u.differentiate('y',out=uy)\n",
    "v.differentiate('y',out=vy)\n",
    "S.differentiate('y',out=Sy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we set integration parameters and the CFL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "solver.stop_sim_time = 2.01\n",
    "solver.stop_wall_time = np.inf\n",
    "solver.stop_iteration = np.inf\n",
    "\n",
    "initial_dt = 0.2*Lx/nx\n",
    "cfl = flow_tools.CFL(solver,initial_dt,safety=0.8)\n",
    "cfl.add_velocities(('u','v'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "We have a sophisticated analysis framework in which the user specifies analysis tasks as strings.  Users can output full data cubes, slices, volume averages, and more.  Here we will only output a few 2D slices, and a 1D profile of the horizontally averaged concentration field.  Data is output in the hdf5 file format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis = solver.evaluator.add_file_handler('analysis_tasks', sim_dt=0.1, max_writes=50)\n",
    "analysis.add_task('S')\n",
    "analysis.add_task('u')\n",
    "solver.evaluator.vars['Lx'] = Lx\n",
    "analysis.add_task(\"integ(S,'x')/Lx\", name='S profile')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Main Loop\n",
    "\n",
    "We now have everything set up for our simulation.  In Dedalus, the user writes their own main loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2020-10-30 21:58:16,584 __main__ 0/1 INFO :: Run time: 152.722414\n",
      "2020-10-30 21:58:16,586 __main__ 0/1 INFO :: Iterations: 269\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmYAAAEvCAYAAADvkw2zAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABGP0lEQVR4nO3da4xk6Xkf9v976tb3e/dMz8zO3pfe5UW8ibYjI5Zj0yCZOGvD/kDJYWTDwEJGGMhAkIiJAeVDvjBfDDuAHGJBE5Zhx/wi2SIEKrIs20hshRIvEkWRK+59dmd3Znqmp+9d3XV78+F5nlN13jqn61RXzfSZ7v8PWJ4+l7qdri6eef/1PK/z3oOIiIiIzl501k+AiIiIiAQvzIiIiIgKghdmRERERAXBCzMiIiKiguCFGREREVFB8MKMiIiIqCDKZ/0ETjLhSn622E+RiIiILrjHP/FhAMB3v/vde9771VHuq9BXPbMo469j/ayfBhEREVGmr3znOwAA59yNUe+LUSYRERFRQfDCjIiIiKggeGFGREREVBC8MCMiIiIqCF6YERERERUEL8yIiIiICoIXZkREREQFwQszIiIiooLghRkRERFRQfDCjIiIiKggeGFGREREVBC8MCMiIiIqCF6YERERERUEL8yIiIiICoIXZkREREQFwQszIiIiooLghRkRERFRQfDCjIiIiKggeGFGREREVBC8MCMiIiIqCF6YERERERXEWC7MnHOfcc792Dn3unPuSycc95POubZz7m+M43GJiIiIzpORL8yccyUAvwzgswBeAPAzzrkXMo773wH81qiPSURERHQejWPE7FMAXvfev+m9bwD4OoAXU4777wH8KoCNMTwmERER0bkzjguzqwDe7Vm/qdtizrmrAP4agK+M4fGIiIiIzqXyGO7DpWzzwfo/BPCL3vu2c2mH99yZcy8BeAkAZlAaw9MjIiIiejSM48LsJoDHetavAXg/OOaTAL6uF2UrAD7nnGt57/91eGfe+5cBvAwAq64WXuARERERnVvjuDD7NoBnnXNPAngPwOcB/GzvAd77J+1n59w/BfAbaRdlRERERBfZyBdm3vuWc+6LkGrLEoCvee9/6Jz7ed3P75URERER5TCOETN4778J4JvBttQLMu/93xrHYxIRERGdN+z8T0RERFQQvDAjIiIiKghemBEREREVBC/MiIiIiAqCF2ZEREREBcELMyIiIqKC4IUZERERUUHwwoyIiIioIHhhRkRERFQQvDAjIiIiKghemBEREREVBC/MiIiIiAqCF2ZEREREBcELMyIiIqKC4IUZERERUUHwwoyIiIioIHhhRkRERFQQvDAjIiIiKghemBEREREVBC/MiIiIiAqCF2ZEREREBcELMyIiIqKC4IUZERERUUGUz/oJEBFddCV3do/d9mf32ETUjyNmRERERAXBETMiolMadaSr5M5wqCx+DrJs+/EPnXE0jmh4HDEjIiIiKghemBEREREVBKNMIqIeeePJk2LIrPvIus2wkWja/QwbRfbHjC7Yn31/+SPK0bJMRqF0EXHEjIiIiKggOGJGRBfKoNGpQaNavfvTtqVtj5BcD59D1vZBz+kk4YhXOPqUtT/c3kk5Pu992Shc1ujb4BGx4YfMOMpGjzqOmBEREREVBC/MiIiIiAqCUSYRnUt5Y8FBcWQ1con99q/ZauT69tl6NUKwnlxWIjnA6YOUqqXEeqRLp8dFQ1QHdDTL851Ocj1eJre3W7JucWOjo9s1Euyu2/7uz2H8Gd7WtqfFoWnHhdvTotBxx5+MPqloOGJGREREVBAcMSOiR9qgkbGsL9p3R7tke6Vv1EuWk6VIly6xnC6XUJmpyLHTtqzKfWVsL03IemV6AgBQnqjpdtkfVeUjuVQpJ9Zt5MxEpVL8c6fdBtAdIYtHxJot3a8jYkcNWW/I9tbRceK4lu5v1Zu6rtvrsmz2rMf7dNludJLH6nOpt5MjaeEya1TOtjd9+vasbXnW++UbMuPIGj0sHDEjIiIiKghemBEREREVBKNMInpkpMWWYWSZtW7R5ERGRDlTlvW5isSEtXmJFycWJHacXJTl1MqUbF+excTyvOzTZXVhRh5zfhkAEM0uynJmQZbTswAAX52WZUWiTF+q6pOt6Lp+NEey9FHwUe16/k3tJTZ0ukRH4kS0Zel03bWbul0iS9c81u0acR7syt0dHcp6/UDWD3d1/x4AoLW/j8ae7GvsyrGtgyNZj7fXZXnQ0KU8dnO/mdhuUWjDtjcllrVoMy0KzYpDswoNwv2njz7zZ5mMPWkUHDEjIiIiKgiOmBFRYWV9cb93f1bLChsZs5EwW85XdN1GwJZlBGz6kixn1mX0a+qyjHrNXF2V+720Lo+zelWWy1fQnpIRMT8pt2lVZcRsrynjNoe6PNIv4B+3vG6XkaHDelu363G6bOpoUMfLSFKnIyNMaSMx4Tkq60hgSUfVKiV5nbWSjvTpeajpckJ7e1RndL2cPI923KRur/kWJhoyMuaOZRkdyWgadNnZ2ZTl3pY8b11v7cro2/H2ni73E+s2Ane8e6zLRmK9sd/sjrIdyrmpt/MVGoQjbFmFBdktPdBn1NE1jqxRGo6YERERERUEL8yIiIiICoJRJhGduby9yML1yVIU9x8Lo8ol7aY/tTQJAJi9IjHj3DX5Av7s9Uuy/oRElNWrjwMAKteeBgD4BdnenpXjDiBf0N85bsfLe3sSp93fkMht9+iOLI8lbtuxuK0ht9nT7Q2NLOvaU6ytMZtFmrYeLk24DgClyKWu27Kq56ccrFfLcp4mtehhSs/bpC5nJuT/JmZqspyqlDBj+6pyLucnJNKdnddjVrTXm/4upnRZa0skOXkkkWZ0tCNP9mAbQDcCbW9tJNaP78tx9c2dvvjzaEsKDdJiTwBoaqHBkUagYeQ5uNdat6jgrAoKGHleLBwxIyIiIioIXpgRERERFcRYokzn3GcA/CMAJQBf9d5/Odj/NwH8oq7uA/i73vvvj+OxiejRM2x0mdWDbKYcYUmrChe039jsukSWC49LpeT801dk/bnHAACV68/JfV+VZXvxGgBgP5KqxU2tlLx7KBHYrbelUvDeocRpd3alZ9f2YRM7On3Rth5rkeWxLts2HVIrOXF4x6oubbu3iceTmZWtd1KiS+9PzrecSz/JUWQTpCeXkUuulzTqtPsp63kulSLUgrhzVqeasrhzVuPP+SnZvjglMfCsHr84WdOlVLjOLUiMPL2qv1d9rOmSvMZpjTxn6zuIDu4DANpbd3W5kVgebe7oUuLSI41Bj7ekgrS+Jb+/buWnRp9Br7WDlvwOe6POQRWfZ9lDjXHn+THyiJlzrgTglwF8FsALAH7GOfdCcNhbAP689/4jAP43AC+P+rhERERE5804Rsw+BeB17/2bAOCc+zqAFwH8yA7w3v9uz/HfAnBtDI9LRI+YrL5kg0bIsr7YP782jXkdGVt8RvqNLT3/BABg8pk/BQAoP/khAEB76ToAYBtSDHD3UEZEbt6WEZT39u4BAG7tyPqGjqTcP5DldvBF/lazjZZO3m0jYzbCFa6HI18+HiFDcnvWiJlP3z5oWx4uSh9ZyxpJc5HrGXWTY0s6ghmVo8R6WQsKbLTN1q2gYEFH0uYnK7ouy6UZ2b6gI3GLun9lag0LczLKNrOio2v6npiO5Pcydyi90xZ0Cf29xiNsm7cAAI1NKSw43JDjbIQtLCrojrA14tE1G1U7inuoDddLbZSRNRYSnH/j+I7ZVQDv9qzf1G1Z/g6A3xzD4xIRERGdK+MYMUv751bqNblz7i9ALsz+XOadOfcSgJcAYAalMTw9IiIiokfDOC7MbgJ4rGf9GoD3w4Occx8B8FUAn/Xeb2bdmff+Zeh30FZdjYOuRI+w00aX9uV++2L/qn2p/LJM/r341AIAYPn5a1j64JMAgNpzHwMA+KsSYR7NyZf+39mX/lU3NiSSekd7Y1lkeWtbln2Rpfa96mhM1Womv8jvOz6OGOMv8WdElSaMLAdFlWE8Gd5/2jHhsYP4Tjux7qKT/0Ec9USfYQxqhQJhgUEYdd7V7eVKcnscfer2qvVOi6PPCpampXBgWeNOiz2XNPZcnpIikJUpKSiYW3sKADD7mE5Sb4ULR9sAgGmNPK2ooLV5GwDQ1mXrvkSg9Y0tHN6V28QFBilxpyzzRZ5h9GlRZlrkyd5pF8c4osxvA3jWOfekc64K4PMAvtF7gHPuOoBfA/AF7/2rY3hMIiIionNn5BEz733LOfdFAL8FaZfxNe/9D51zP6/7vwLglwAsA/jH+i+qlvf+k6M+NhEREdF54gb1wjlLq67m/zrWz/ppEFFOWdGlsciyd0oloBthWrXlaq1bdQkAS8/KlD/Lz0tB9/JPfAAAUHv+E2hflp93qksAgFsaH71xX/qPvbMjcdN792W5sWeRpcRO+xpZtpraeyzoOdYbXQLJnmNZfcYGVV32RZQ+fX9adGn7+x4ziCQHrYd8++T9J3GlZPwZxqFZ6xZ5RhZ12rpGm1HcQ02PL7n4Z4s9q/peqdSSvdOWZzTynE5Gnra+ohWha9NWAapxuVX81jRObUr/s+hgE5HGnq277wFIjzsBDIw8j3Q5qHfaST3ULNJsZlR2mkHbs4THpx8z8JAL4yv+bQCAc+67ow48sfM/ERERUUHwwoyIiIioIMYyJRMRXWx5qy8nS1Z92Z1SCehGl8uLUlFnVZcrz1+W/R+T6ZMmP/iTAAB//cMAgO2JFby3JzHQa7e2AQA3tiU2uqnR5d297hRKQLJBLNAbVcpztdgwq6rR4rVOzz9rSxrBdSNIO1jWI+0qFBRpDt0UtjfitGgyK7IMo8msKHNgxDlgf5pBlZ0Do06NRm27LUvlcl+lp4srPOWYu/qeeruSrPS0qLM2qZHllDW11ahTI821OYlAV3R9WfdfmlnB0qx8tWZ2Wd5/FnfWGtKUdlKjzpV9aWobRp6NuzJtVNjU1iLPw3vynj0Kpo1qHDTR0PfvfitZ4ak9jnvWk9FmN+rM28xWNtjfbFpc2T2mf18aRp7D4YgZERERUUFwxIyITiXtX8vhSFk16Es2XUpOrXTZ+lNdnwMArL6wAgBY+cjTAIC5n5DeZKUPyEhZffEJAMC7uzJ68Pq7u/GX+9+5J1/2D7/cH4+QtdKnTwpZL66SvggfpfcScx0P64Edj2RFPft6trd1rCxr5CxrSMEHRQHdIoJ234hZ1gjZaYsCTjNSNqyskbWTtoejaFmja91RtuQ0UYN6pllxQVXfmzM9xQThqNrarI6uxYUEMi3Y4oKMrM2tfRQAMF8LeqcdSC89tytFA60NG1mT6aKON2TEzUbW6ne3UN/cl59tNG3HJmFP9k6z9/s4p4k6bc80jqydDkfMiIiIiAqCF2ZEREREBcEok4iGkvZF/0FTLM1rbBRPsTQvX/K3/mSrH7oKALj0ky8AAGof/rMAgPbVDwEAbrfl+LduS1z52n3pLXVzqx5PrbSj0eWe9iVraHTZCqZLsqjSlWzaoLxTF/VHmlkFA2E/MvuCuhUaOP0ncVueKqJOsqDAB8/VB9PoJKLMdhBpnjLC7JxhpJklLdKMcvZIy+ydVqkm18tVvd9kTzWLNC0CvVkrdWNOjT3T4k6gG3GuWuRp2y0KnZb3+/zqdVle+xQAYLYsv+cZjTrndJoov3W7p5BA4s7D23JMfWNbllpAUN+Uv42snmnHexrxD5geqnfZX0DgEusmLCw4qYAgebv8kedFiDs5YkZERERUELwwIyIiIioIRplElMtJvcqy+pRZ9aX1KVu9PAOgG2GufewpWf+kVF9Wnv/TAIDjNZlm6R2tOHtdo8u3t3WaJa3A3DxoYP9IIppDrUZrB5FjuZzv35/Wi2wQu//ex4l7oVnUE/RGizQ2KpWSlaHOadxqd2S3s+eSo89ZVoSZFWl2MraH99e3fYQpm8apUxq+krNXNCjqLCejTqv6jMpVlMryf5kWd5aCCs+b+j63yLM2IRFn2DttbU6ieZse6vKCrK9adefMJADg0syTAIC5689g4Tm5z4m2VCEvaty5vCe90bJ6ph3c0sjzrlZ43pOeaxZxxlWeJ1R3jqOys3cdwXF5eqaZQXHneYg6OWJGREREVBC8MCMiIiIqCEaZRJQqKzJIayJrVZjdBrISu1yekOXCukSY1kD20idliqW5j31S7usDUpW2MyPVaje3JE75k3tafbmTnGZpp27TK7XQ0hgljCKrQYQZ7i/b1D7Bsu/1aizb7iTbwrY7Po4z7TlYJWhYEWrRZUvnz7Fo00XJeLCFpLi60/dHm4MaytoyK7ocduqmvPsfOD1Jg6Z9MuFx9uxdKV91Z2/0mdXE1io9LQYNI0+LNi3yvBE0ta1O2nRRVt0p97M6KxHn2lwNlzT+tMrOSzNLcuzCJQDA/KWPytKa2eq0UFbh6e9LNadVdVr0GVZ3HsaRp3xdoL51FMebFneOb3qowdWdD6KyM8/9nSWOmBEREREVBC/MiIiIiAqCUSYRJQwTYdoyK8JcvCZzYK5phLlmEeYnpPrSPf0JAMDOtMwt+K5GJa9uSoxyY0uWdzVK2T+25rHdOM0iyaxo0qLIcHu1lL49vJ8s7Y7vizAtTrGKNnueVjG6r81vW01Zj+dv1HMbabTZ0MeIgsa11gS1N1TNiijDCHPYOTWz7j/9mL7ZPx8Yp7/PYSPVvuizNWB/yvZwXs7wmLDis6+JbdDctlyVxrNh5HmnYo1sdTlZiSs7pyaTlZ1hM9vLs8nIc332CQDA/HWZf3b+WXmsyY5UY2ZWd+ryeONuf2VnMG/nUVzZqZHnjryDjzTStMjTos2sKk/7++mt7hxHZWfvcabIzWw5YkZERERUELwwIyIiIioIRplEBOB0ESYglZhZVZgrH5DKsdWPPQMAmP3IxwEA0VPSUHZXI8z39qTK6437El3e2pNo5P6+RCJh89huPNnzPO15aTVmPG+nrtd0Ge8Pjg+3x0uXfmLa3vc1m7VIM27Oqcs9jWCtGa7N57mn1W02d6bNqdkJGst2qzNt7s3+2C1vhJl/Ls3sePIsKzNPG2EOjlvtd5E+XuGi0sD4s2+9lBFt6rIxIPLsrfK0ik5bvmfzdmrEaRWeE9rEdkErO8PI88qCNK/Nqu5c1OrOOZ3Xdnr/Lma1wtMqOy3ubG68D6Bb2Xlo83baMog6bb15oFGn/h10qzv7I87BTWxdsI6EMOrsbkeq3sgzb4XnuHHEjIiIiKggCj1i5pB9xVrkHiREj5LTjpTZtEsz5QhL+q/r+bVpAMDSs/qv8A89AQCY/eCH5L6e+SgAoD5/DQBwe1dGKd7ZkX9Jb+i/pHd0JOm4lRzlCHuTAf1f4reRsUkbUQhGxCZ1xKGix1f09rZetvvT110pZf/71UbKmtbPTL/sbF96tudvfdfsdW3r+pQ+x00dGTyKZLuPR8xkRKHT0l5qLet/Vur24Qqbn6lBI2XdZSf1dln3N06d4D6jnL3JhjHoeYejXNnHN1NGKpvBfUXBer5eaVk91WxZKlcHjqqVqzJSZiNqVkDwphYQWI80KyKY1hEz65m2Pi8jaTY91OUZGWFbm5nEpRkpHFh8Uop35p/X+2zI9E6T+3cBACvbMh1UayNZQHB4R0bcDm/dBwDUN3d0KX0KD+9Jf8K0fmkPuldadz1fEUHafY0bR8yIiIiICoIXZkREREQF8QhEmenDkMN8KY+xJ1G/vBFmd12jP13OlK13WYQFjT/mrs0CABafviz7nnsSAFB96oMAgObidQDAvUPJ327vS3SxeSjRRbdPWTJei7+Ij+QX82vlqD+qDCLMiZJ9+b8U3wYAJqwowKJMfcGVKH39hEQTzbZFmhpl6vM/bMpyUWOknSl5fVv6ejesX5U+t43gy/8d6+uk929TO3XaHYSB26Av9WdFmIOKAEwYOw4j7J2WpT3guDDyG4uM15UWqw4bi/ZHnWFkmS/6bKE/9kyLO3vXw8izVJWoMi4i0L+T97Ro4JUg8uwtIrCYc13/zi/NJnumretUaouX5e974XG5r2ntyLegUedSRq+0o1t3AAAH8fRQ0i/tcPOgr3DA4s7mvpzbA+0VaIUDWT3STuqVJsvsIoL+649800SdFkfMiIiIiAqCF2ZEREREBVHoKBPojVvyZZdp1RLnYbZ5onEZNsIMqzEnS1aNKeuL1TJmtBpz9to8AGDuSelPVrkuVVzt+SsAgB2tLtw51ipF7WNk0V876N9VC6ow+6ZVKkdxdGnVlrW+CFOWU8F+W7dIM4ws46mZ9HzYtElpbc3sY8c+fxrt5HK/Ia9zpybLWXvO1eRzMnYevOYtHT0/rYYc14pc5hRC8XPKrMbMF2FmRZd5Y8kTn1vOWLSvYnIMj51XGKvmilEHxKJZr7v7OrOjz0GxZyucDqqUM+oMIk+bJiqu7qyV+yo7J6ZlOTud7JW2Pq8VnQsWccp9rc9m9EqTu8GMRp3z2i+tffemvKa776FxR3qn9U0LdU8qQg83paIznBaqoVFn3FOwnd0rLW3ZwfCVnePCETMiIiKiguCFGREREVFBFDrKjFw3PslqBtdvcG6Zddthp19g9EmPgkHv6+wqzOTtJ6KwGlPijamVSUxfmpJ9V1fl2GuPAQCiFWkk25xeBgDsH0qcsHcs8cKxVRkGf5NpjWR7t/cuLbLMiirDdVtaVFmNqy9dYntFn0IpiDRtvZclsDZlkhZjxlHmlN7ZrFbAxc8haI7bndpJzk9DY5jmscVPsu4ypokCTqrO7Jy4P4wuT4oNH9aUTGc59dM4DYxFh6gMzRN7ynqU2N4KjhsUefZGnYMqO9/R97VND1WbkIxyetaiTok0V7WK89qSNrOdsepO+SrEpWmJPBeeex4AsPiREiYPJN6csWmhNqWSM++0UIf3ZJo3a1qbNS1Ut6qzG3n2V3RCl/p3/oAazXLEjIiIiKggeGFGREREVBDFjjLhMK1D/YOawSHej2A9bajx9BWevdjkloosb4SZdTv7GkElqM6c1AOmNbaYWJzA9Jo0lp1el8iytCxVmZ2pRQDAUUduc9SSP4Q4wuykR5hhXNhbhQl0KygnylFcjTnRF1km121OzTC6tHVbdqNMWUZe48OOhkEaMzqfbIILAN5pk9uSRDntiiwnW3Jf9VZyXs5K8EuyxrqHGmHuaczS1Oi3dCTxVakcDazKjJ9TzurLrOgyT5x4XiJHIK1JrBjna7T7GvQ7TGu4O6g6NIw/sx6jG2kmx2f6qjt75mW1eDOs8LSIM6zsrGjE+b5FnVrdOaHVnL1NbAH0NbK9tjCJNT326pw0q1548hkAwGIwX+fUnjSp9felitOa11rkaVHnwW2Zr/PI5uvUqNMiTqvubBw00dC5bbPm62x0GGUSERERnWu8MCMiIiIqiEJHmdXIYVWHQLPmu+ptBte7vRsdur5tYcQzSoXnybfvGrbis/8xRrs9XRzDRphZc2KGDWa7jWU1jtOIoTZXQ21Boszq4gIAIJqVZVsjDqtOjP8W9e82CmJDaOwY/k2F81f2No+1OTDD6LIbYQaRZUaEaVFnHF22JMZwbaneQkuWcaTZaffFmRZlwmIgjXymK1K1Wq1JxGNNbO1XYXNtHmtksj+nDWnr8ph7WkFm1XAushq7rjCSHNRQNivCzIrszlNceZKH8TrdgIazJxkUg+ZtjJtW8Zl2vy4qZTa3jZvVVoJoc0DUGVZ33tb/n3/Tqjp7Gtku6LycYRPbxxblPta0svPqrFSCL119AgCw+LQ2ce5INLmoUefSjixbG8n5Og91vs64mnNjC/XNA9kWxJxW4WlR57hxxIyIiIioIHhhRkRERFQQxY4yayU8tSYRQFPnvTpqJ6sjBs171fbdWCQr7uyv5JTluCPPk+/jZKNEoYxBz78874/BEWZyu/2rrTtXpqxbjFidkbihNldFZU7+TqNpiTRdTdZ9SaKLsHipFESYFkdGnWS1ZjfqTEaZtVK3KjOMLicqyerLbnQZbA+WUUvm2HNtjSeaR4n1ONLUKNO1W32NQbtT+8pj+bKcI1+R+67WpJHmXG1Otut5aU3JedrTasz7Wq02o406uxFm/g+CrGrMvBHmMDHbRYk581bCDvIwz1f4WPYa0io+gfToc1Clp2vo30pkc7n2V3b23nfcvDaMOq3qs6ZRZ7WCOxpzvqVLq+S0ZrazQdR5TSPOq4vJJrbrszJP58qqNL1efOxTAIAZyN/1gs7XubS3Ia/hzrtobcicnccb0tz2QCs76xsyX+ehRp3jxhEzIiIiooLghRkRERFRQRQ6ypyYn8TTf1maydU392XZVx0hEUFDo866RgHdRnC+L+60aC+s8Myu7MSJ63kiz+6xw2WSp40+e41aERpiNHr2hvmd5o8w06sxw6rMkmaa5Qn5+ChPVlCekAgCGkX4KP2jxe67FldVWhWm7rd0RZOSsGqzYrdLiTJr5aB567ARpi7j6ss42rR1XVrU2WnBt+znjLkOqxKv2HFeY1A7OzNViX6PtBpzXuMam1PTmudaY91Iz1/U8zsN57580BHmRYkt0zys154nMn1QzyUt+swbew6MPIOmtn1RZ1DlWSpXE3N2Ar0xZ7KJ7bu6/GNrfK1fDwjn67Qmtld1vs51nb9zfUbm6VxbkMhz8dJHMV+Wv0ubp3Ne487WxrsAgPadd1LPy6g4YkZERERUELwwIyIiIiqIsUSZzrnPAPhHkADiq977Lwf7ne7/HIBDAH/Le/+9QfdbWV7F41/4WQDd+a4ad2QerINbWh1xV6oj6vdkviyLOOs9jeCa2pzxSOeds2gza96rMOLsVm8OH3kOX+EZ3u50OeQ4ItAs445GT3LRY9NRznVWhDnosaJ4PWw4q/urVkmly0oZrhQ0oNTGq9aMtVKSuMDiQ4sfm1qFaffd7Mj9WAPasHozLcq0CDOrkawWafYsNQ6Mqy01jgwjzJZWmgWNZn1T9vvjOhBHlBr16OuBRTXHeh+1ieR50ai3WpmMXwcAzNaSFaZTdo6DasyO95kR5KDGsuHxWeuDttODM85zHsaJw96377Rzzx/qs5rWjhB5xvGmNa3NmLczjDqt0rOqEed7uvyRRZ3BfJ2XdH7Oyxp1XluajCs6r87JVw4uzchcwAsrH5TlRx/M2NbI9+qcKwH4ZQCfBfACgJ9xzr0QHPZZAM/qfy8B+D9HfVwiIiKi82Ycl3ufAvC69/5N730DwNcBvBgc8yKAf+bFtwAsOOfWx/DYREREROfGOKLMqwDe7Vm/CeBP5zjmKoBbJ91xPaph47lPAwCWfkLnvdLqCKuS8JsScba0OqJ5T+a7OtRGcIcb2/HcV5kVnQcSR9i8V8M0r02uW9TpEuu9x+Zdt2gnbyTZf/vhM7AHGX+e1sOMTR91Wb/z8BwOqsYM91tjWdsenfRLsWjD4kBt0lrT6sMJvW1Dc8WOl4+gSiR/cxZtWkPaKHiO3QazVq3o4miyHC5L6dttLkxYhNlOVlt2qzB1XSNNr000/XF33bcs3gyqMy2KqVT0PuS4KN5e0/NTl9dRmtHzkKy+zOI73SizY8+hb47M9GrMQRHmKDHaRY89x9WAdpxO+zvpfS157yNrHs9xRJ4usia2w1d2Aj1Rp0acFnlWarLd5ut8pSfqtJhzKajovKYVnVd1HQC+8PFrqa/pNMYxYpb2CRL+P3yeY+RA515yzn3HOfedrc3NkZ8cERER0aNiHBdmNwE81rN+DcD7pzgGAOC9f9l7/0nv/ScXl5fH8PSIiIiIHg3jiDK/DeBZ59yTAN4D8HkAPxsc8w0AX3TOfR0Sc+5470+MMQHg7n4D//R7ElWuz0vFxFVtBndl7kkAwOKT0oB24XltMnkkVZpTe9IIzt+/FVd02vJQKzoP724DAI42dwAA9XuHssyIOrPm6+xGnAjW/QmVnf1xZ+9xedfNsNFn2n2dtgI0/b6LF4ueRyf9zrIizKzjomA9q6rTIgTTaXfQPtI47UjmjuscSJW0m5GortKR/ZMa7bWDqNIqK49bVvOcFAVRpj23cgl9UaZVZXYrOXW/Pf9WGFnaUqOSZrIqM44w6/raeiNNixEzGs16i1GsQtIqy3TOTLt9ST+Jw1+RVac29Ly09TOn3erE0WRnyGrM+Llxbsyxe5jn5UHHpqd5Lfachn3vZd0uLfLMW+GZFXlmVXeWU+brrEzI9YbFnG8EzWsnpivx444zyhz5wsx733LOfRHAb0HaZXzNe/9D59zP6/6vAPgmpFXG65B2GX971MclIiIiOm/G0sfMe/9NyMVX77av9PzsAfx343gsIiIiovOq0HNl7h828I3/dANAd8hwblarIxa0OmJRllc06rysUee1uccBAIvXn8bCc1rR2ZR4ZULnu1revg0AaG1IxNm2qPOOVHwe3roPAKhb1LkpUcbxrsQPFnlaA1ubr7PRtPk6fRxrWoWnRXxhZefokac25OzJQgY3Zz1d3Jin6es4Y9GTXJTIdJjzOSjCzKrGDI/PijQtKutorNZpttCqy99AWyPM0v42AKA8L+teG6lOT+v3RrXU0x7jKJLfo8WQnazIPqgQLbmeKDNsLGuRZhxhBo1jm0fBdqvO1HgyjDA1pk1WZSab01q8GMctFitqnOJ17sxIG++G7P3c0nN7qHP/NnTZ0s+WTruDdkNiYqvKjKszB1RjnjbCZHxZLOP+fYwjGh32OQ2KPnu3DxN3AtmRZ1Z1ZzOlqjNsYhtWdNp8nePGKZmIiIiICoIXZkREREQFUegos3XcwK03JW4M5716e1KiTYs4rRHcbBB1rs9P4LGlKfnZ9s1cAQAsX74OAFh6QqPOjgxxLmjUubi7Ic/jtjSvbW9KIWn9ljSxrYdVnZv7srSqzp3jOPa0uLN5bDFnsqKzG3UisX1Q9GnyNLXtbu+PPU86vt/o8eG45sB8WJFp0aX1fB0cYZ58/CBef4mtoxaahxrr70plc2lHKp+jKWksG5Xk79M7+bfgzOQ8AKCsf8/VljWYlfu292iYVLvwNUT9VZmVoJFsWGU5aD2uLNUI0x/Ja/LHEh3GEWfjKG4sa41mYxavWIxY1tcfRjYWo+jrPtIf9jS63Knr1yN0nt/mkUaajeNuhNm0CDNZjTlqhMno8mIZ5fd92hh0HO+xUSo8gZ7qzJTqzjDujCs6gya248YRMyIiIqKC4IUZERERUUEUOsrstJo42JAYcdC8V+VJmWvOos53ajbfVSVuBjc9K8s1nffq2qJEnFd13qv1GY06Z1cAAKtLEnkuXf0EAGDGSawwp/N0LmrkafN0tjcldj2+LcvDjS3U70rD2+MtiUXC+TobVtGZ0cS2G3VqY8kg6jypurM/7kTi2FBWpWfIYqTTxJHdGHX42w7/WA/+MR6WPOdrUAPZrOrL/sayJ99fpyfCBIBWvYXmrr6/NdYvT0v1odMqxLJGeSWrNtaqxMma/N3WtGqzpc+mpe/rrOrMeA7NyHVfh0WXYQPZoPoynhvTIkydrzKMLuP1RjLijKsyWw10mvI6fDvZGNdp+WgcjgRzCPpIPp98Wc6PNdbdOZb7u7cvz21Tl8cWaR7bOd9HK6jK7KvGZIRJD8mo75k8UeiojzEo8uwE24G0ys6wovPBNPnliBkRERFRQfDCjIiIiKggCh1let9BY1+iwO48V5XEelgdEc57VapOxs3gKjXZd1Njzh9pRWdVKzynZmT/pQVtVjuvTWw16rwyJ9vXZ6Si7NLsKgBgYfXDstRps6Y14pzbv4f2pszVbs1rm/ekovPwtlSt1Td3AfTP13mk83SGVZ2tulZnZVR39s7bGVZ2ZjWzNd2mtia9uW33+NTNJzZ9HSUGPUnaYz6MuPSsjDJH5mmrMe131tGKwbY1Pz1o4nhXYrXytFQmlyf077GS/Igpa8xY0jix05B40FWn9Xj5G6uUtNpJqzj7JpE0rTacNWvVZbdRrDZ/DaNLayCr0WQ7jCyD9TDStPiy02yh00hvFBtV5XWXatqA0j6vNNr1+vnUrsjXKXYP5TltaHR5c0vO5/aOfvVBv+rQqNf1JdXjasyOntN20OQ2xAiTimqY996g+HCckafrJB8rLe58EDhiRkRERFQQvDAjIiIiKohiR5mdDtrNZONGFzRyjIJI05aNnoZwYcyZNe9VZUJihvc16qxNamWnVnVOatS5pFWd60HUeVXX12clllmfWcLCM88DABY/JM+rVpf5N6f2JO702xJtWtTZ0qVVc9Y3tgF0o86jbYk6e5vYAkDTog5dtuqtvsrOrErOwZFn2IhWttschMl6tPSYMtwWxmmjz3lZjNwy7+sYV2PcYRrLhvvDf5VlzZ0ZsqrMdkN+8616K37/HW/rXJkaYVp14qQ2P7U5JSONEaOZBTl+eg5At0rRVywC1OpFl/5vSOc7gNf7jqNMrcq0iFOrMS2KHBRdxkutvmwfy+3bRzofaEObynY6mdWYFmWGEWY0LQ13W9pgd0e/knBbI8y37svf94178pwO9SsN9X37O5fPgWZ9H217ngMayuaZh5DoUXHa+ThHvZ9x3TYPjpgRERERFQQvzIiIiIgKotBRJuDj5okmb5VEd/6ro76Y0+a96lZ6WgXZoKhTKqhuadT5+mRynk6bt3OuZ77Oa4sady6GTWwfBwAsXnsaALDwjM7X2ZIIo2bzde5JM9s46ty4CQBobEokerihkWdKI9tjjUHi5rVDVnb2R53Q5cnzdvZGn4Njz/R5O7v7Uzf36Y3rRo9Fh9P7HE8bUY5aQZr2uIMizHDOzEHPIZ6/0qLMprx/mketuIrY6Z3Y35ixSsYJjQMrBxJ5RrMLstQo0+nfmEV/Lv4aQiXzefkgJrX5K+PGqxptxo1h+yLL5ByY9lzD6LK3GhMAOj0xZmQNZUs2l55+xtQmEq+zMyPNqxsTiwCAu9vyGG9qhPnaHalq3bovz+nAqrP3tgEAzaN9fU71+PW1m8kGs4wwibryvt/TIs+z+lvhiBkRERFRQfDCjIiIiKggCh5lZs/zFg479s2DpZFnJ+XY/vmv0is7w2a2FnVa89pB83W+PVmJKzutonNWY841bVZrEedVXV6Oo85LAICV5WsAgMVgvs7pA4ky5zTybN/TRrabtwBII9u6xpyHd7cBAMdbEh8dbUtMYlFnWNnZ7JkDEQAaGlnlr+7sjROTDWVPij3T9ueN15K3GS0WHdZJz7H/9TyY6tE81ZnZEWZ6FWfW6+pYI2JtMNuqt1CqyL0fl7SZq0Z7FjO2m/aekv3VuQNdSoPlyoxUMlv0Z1Em9G/MWZR5UmPHeM7IINq0KNMqvK0hq1ZbxtGkNottHQXbNbJMayZrrxMWZVbkeZan9GsQi2vytFfk77g1dxkAcPtA7utPtPryB+/JeXjrliz3NMo83JH1xqFUY7bqEmV2Ws3MOTEZYRINb5xNbkfFETMiIiKigij2iJn3uUfI+m7ac7us0TQzuIAgOZLWjHZ1vdsrDeifFqpUrvaNptm0UDd0VO0H01Y4INutkGDwtFDSD2l1ehkAsLj8QQDAwoT2Szu8H/dKW7JeaZu3E8v6Hdl/ZNNC3Zd/lfcWEABAQ4sGbHooG0lrBSNrzU63iCBrVG1QAcGgEbRw1KsST/HUPxwWFhqcZvRtVA9uhCz7frP6k3XXk9uzRs5C4e+q1FMEYO8FY0UAPh5t0mmcdDSqeajvrV350ntlWvuf6VRO5UkZOY6Cfmhh4U6v7ohRJ/nYwYhYR/t9hV/mb+v+8HZhr7Lua4ziL/vb86zMSfFCaVlGxsrrTwAAWkvXAQC3j+V5/+iujHz9wbvbAIA/uiGj2/fvyN/e/rYsG3syMm4jZda7rN1qxFMycaSM6OFiHzMiIiKiC4IXZkREREQFUewos8dpe5Gk3a4v2tTbnFRAkHbfvb3Sevf39kvL6pUWFhBY1FmuSoTzvvZIC6eFsl5pCz290oCePmnzFnXW+nqlzT9t0+PI850/2JT9urTpoCzqbN6VCDSeDkojz7quN/aS00JZ1Nk4aPTFnDZ9j00TZXHYsJGn6d+eHWkO+tK/3SYaMnZMD7gejDz/gsqKIgd/yT890syKfu18Wj+zTqONVnCw9fWyY8LeZ+UDee+0puX92tiT97NN5WSRpkWYJZ3iKI4yS9lnxKJHizTty/tZEWW8PYhAbYojY3/v9tjlahllfZ61BflqQXn1KgCgcv05eX0rTwEAbrfldf7orkSUv6/R5e+/IX97929LVLl3X5ZHW/I32NApmFo6hVVbCxo6zQYjTKJziiNmRERERAXBCzMiIiKignhkosy88gzf563wDCtC+47TQrROGHFalNnTLy2s7GyM2CvNpoV6IyPynJ6txjGnxZtXl5LTQq1NL8lyTnqmza9+GACwUJXr9QmNOKd16TOqO4/vSeVYPY48d3C8LZFM3CvNYs59rew8Sq/stMizt8Kzd5ld1dndnhV7dtcRyI5B09jtH2wnm3xOqjQdFF1mbe+fyin9/q2fmWv7+PfmIv29QSO3djJOrPT0PgO6v/9SRfp2RdWSrpd1PYwwk3EiAES6z6LIUBhN+uC4sOrS1uNK0FLyse25VaYnMLkmUytVLj8GAKg+JdXRrdVnAAC3WvK3932daun/e0v+Vr71mky1du99qUYdJsKU19BmhEl0TnHEjIiIiKggeGFGREREVBDnLsrMY9gKz2Ga2SaOb3W3PcgmtkA36izVZFmZmMDNCYk5qxZzBk1sl+a0sjNoYruuTWzX9PjM6k4vscqsTgs1fyiVZu1778dTQ7W2ZJ9ND2URp1V6NnaT00Md6TJuWhtMExU3Ks2IPHung+puC5vW2vLk6aG629Mjzgc1xdNp5JmSKTx20FRMWftDvu3j+M+qLrtafccCQKmpx+vvs1S1iNIl16Pkugkbzp4kjC4Hif/mdN2ax3YrMOXrBNNXV1G5JpFl6bE/BQBorD0LALixL6/rD27J+/w/vi5fB/j+m7Lc0kayu/fk7+JoR/5OmhphWiPZ1nFdX0P2tEuMLonOF46YERERERUEL8yIiIiICuJCRpl5DYoIhmlmO3plp0QfYeTpgurO3sgzK+a09Vta0flqzaJOnbdTo85ZbWK7ZvNzBtWda9Mahc5IdefqjFZ3Ln8Q8zV5XtWjbblPrex0+1KVZlFn3Mx2U6rUwma2x9tStXYcR57dJrZAf+TZbrb7Kjw7fdFmeqWnyW5em35cclvfpocmK2ocFEn2rw83h2avuKFsyaJK26NNXe2cBxGmRZsWZUb6O7TIMtLtLl72VGUOmgRVWSwasma4VhFq6+VpeZ9XZ2X+yymtwJxYl/d55fHn4a98AACwMyXb3rgnUfx33pf38f/7qkSUb74j61tanXlwX/4ejvXvIYww2xnzYDLCJDr/OGJGREREVBC8MCMiIiIqCEaZIximme0olZ2924ep7uydsxPor+w8yFHZCQA3gurOmkagWfN2rs9P4IrGnmszFndeAQAsLl8HAMxf+QQAYLYs0dakNbPVyk7sSbTZuidRZ2d7AwDQ2JTjLPI83tJGtj2Rp8Wd1rw0q6mtRZ0Wq7VbyXk8s5vZJqs9h4k0H0bUmR1phuvjaUALAB2rtiwl1yNYc1c9V3quLYq0aNMiUBdEllHc3DWYizPKH2mG8acdV6ro34xGmFZ1WZ2dBgBMLM8BkOpLACivPynLx58HALRWn8Z7R3Kf378h70ebA/Pb2kj23vsSye/cOwSQ3UC21ZCovpMRYYYYYxKdXxwxIyIiIioIXpgRERERFQSjzAdsXM1sw+1hlJlW3Zm3mW1Y2XnaeTtfr1nUWY6b2U7OyHJNm9muzmqF54Iu42a2GnlqM9v5y08BAOYel+cyE8nrn9LIc+ZwW87LjkSc8fydWxtobsu+MO5s7ElsFDa1bWhFZ1zhGTSztXWL2yyOs3Vrctvb3Hb4+TpPPj6PQdWTw1ZtDtoeuW7sOKgyshttxlv0f22LRpztZITpS8nzYLGktwiz2elua6dHlsbWyxPJhrGVOYkuJzW6nLq8DACoXpH3YuUJjS5X5D25ATn+1Tt1fFerL3/vDXlf3rgp0eX2XXmvHdzXBrK7Up3Zqst7cZgGsmnrRHR+ccSMiIiIqCA4YlYQ4x5Zc1FphJ5pwkbWGlkjakExQW8RQbeAQHpA3dTRtEpGz7Sp6eTImk0TdVlH1lZ1RG1tWu5vfVZ6Ss1e1WminpLnNB21UQlG1bAv62EhQXPLponSwoHtZCFBY1e+sG0jak0tIrARtLhfmvbgatVb8fQ/VlgQT1UU9FQbNII27AjbSQa1+Rrc5ywYeXLJL+gPY9DIWRRv1TIXfQtHpeS/IW1/73YbKStVbWlf6tf33rS8l8Iv90+tywhZZf0JWV5/DgDQWpb1jWgWAPD2toywfv+2FKb8/lv38eq72wC6Uyztb8vyWKdYOt7Tvn02QqZf8vdteWGdVtzoTbZzpIzowuOIGREREVFB8MKMiIiIqCAYZT5ixhFtnLpnmn5ReVARgYuiOObMKiQIe6aVqxJV3qgle6VZ77RhIk+LOy/NLAAA5i5L3Dl7XZ6L9U6b0J5p03XtnXawDUAKCACgo8v2jkShDY06j4Los7kr8VXz4Lh/qqi4kKCTWA8jznYzWVBg+ztBhumD9c4JRQJhHJp3aqUo47i+XmIp0yN1+49lTX+UPbVS+vE2XVKUWHc9vcjK+h4pVWRfRWPvqn6pv7YgRSo2pVJtLdmXLLoi74/WovTYu9mQ98lb9yV+/MGGvA/+4MY2AMTx5fbdQ+xvayFJX3RZ12X+KZbS1ono4uGIGREREVFBjHRh5pxbcs79tnPuNV0uphzzmHPu3zvnXnHO/dA59wujPCYRERHReTVqlPklAL/jvf+yc+5Luv6LwTEtAP+D9/57zrlZAN91zv229/5HIz42nWCUSOQ0PdPC47LizsaAaaKiAdNE5Yk8Le5c1h5qYe+0tVmr8JTlytRVAMDikvStmruskafGZ9WmRJXVuvSsmtOl35XqvPaWxFid3U2097YBAMdbWtm5J5WdTavw3LPYU6v0ggrPZjxdlFbtBVFnJ4g6LdrstH1cERrGnyaMQfPKqr7sjSvj32dwbNZ0SC5Kn3Ipjiir9v5Jrlts2duLrBJElpPL87KukWXpkkSU5Uvy+20tyPRge5UFAMDtAznnP35bfjc/uiO/uz/SyPLmbYmsbVolq7xs7N1H41CnVgr6kzG6JKLTGjXKfBHAr+jPvwLgr4YHeO9vee+/pz/vAXgFwNURH5eIiIjo3Bn1wuyS9/4WIBdgANZOOtg59wSAjwH4vREfl4iIiOjcGRhlOuf+LYDLKbv+/jAP5JybAfCrAP6e9373hONeAvASAKAyPcxD0JBGjVFOjDxzxJ1AvgpPYHDkGa/XJuOmtjeCpra1yWT8aZWeC/G0Ucmoc3XOos6qLud0KQ1J564+CwCYeUKe40y1hOqxvLUnNO50RxKL4UAqP63Cs2PLA9nf3JXbNXaDqFOXzUOp8mvVJSJrH0nTW4s4W/VWHGGGFZ6dIPYMt4d85+TI02LIXmEkGR9bSkaWVjmZ1Qy2VNHlhPxuypNy7q0pbHVuKrFeW1lCaVEjy1UZiHdL6wCA9pws9yO5zcahvCnf2pBz+cbWHQDAD9+Tc//q+7Lc2bTIUs794Y7+brTishnHlvWho0vDCJOIsgy8MPPe/6Wsfc65O865de/9LefcOoCNjOMqkIuyf+G9/7UBj/cygJcBIJpaGaHfOREREdGjZdQv/38DwM8B+LIufz08wDnnAPwTAK947//BiI9HZ+CsCwmA/CNrLiqdOFUUkBxdA4BKTUdlbIRtQu6zNpEsLAiLCpZsmqh4gvYaFnU0bmVqAQCwOCGjOXNrOhn7NXmeUzpyNKEnoHwsI2fTOsIW6dLrl8utqMAfyOhN51D312WErXVYR/NAe2cdaS+1o0Zivd20nmoy2hYXCzR0u42sdZIjaT5jZC3RvywK+4zpeqWs6zY9UjWx3aZJKul26z1WmpYv8kc69VZpXkYpbXSsM7Miy+ll1Mtym/taQHHPRsZuyPm4sbUNAHjtjox0vbkhy737sj8eGdvXSe11ZCz8Yn+npZPd90yrxBEyIhq3Ub9j9mUAn3bOvQbg07oO59wV59w39ZifAvAFAP+Fc+4P9b/Pjfi4REREROfOSCNm3vtNAH8xZfv7AD6nP/9HAMPPdkxERER0wXBKJhqbs4g8w+MtMgOyY8+oXEmu5ygoALpRqEWeZf2iem/0WQ0KDOa1wGB5RuLOZY1DlzQOtejTlvMTS7I+cQkAML2q0ecVWU6WbSn/1ql4jUIbB5jUiM21JJpztt7U6E2/qO6PZdkJ1tGyL7B3EusmjDiNiyLAfi92ru0c1ySqdFVZRpMSO7oJWfqKnFNflS/odyZmZVmT5WFbXud+Ux5771iWdw4kdrx95wjv7W4DAG5on7GbW7Lc3JLXfbgrxx7uy+up78n25oFElc0jiSrDXmQdff0dfW92gi/4hz+nrRMRDYtTMhEREREVBC/MiIiIiAqCUSY9VOPunWbi6LPV3TbOSs/e5UnRp8Wetq9clUivEsSf1Vopsb1Ss+2yPms91rSH2oxuX5hKRqRzWjk6Wy1hRh97qjKjS+2vNqOPrf3HqtpbrKLrtt2KLEtO1/Wboda2zLnkV0W973azsdme2toDzdabut7SZUN3HLUsmpTf5862/OJ2tT/b/SPpvHNf48cNjSPv7knMeK8nnjw6kGOODuW2jboutUdcy6JKjXbbx7rMiCx9O73SMu29x+iSiMaNI2ZEREREBcERMyqkcYxE5C0oaGfsHzTCFh4nPdQqfduA7iibrYcFBn2FBuXkCFukX/q39bKOhpVtcu9KqbvNRuV036Qua1Y4UJX7rtq6ddvXoTE7rhSPpCWXpt0zQ4D93NCRsGNdNrQHWl0nZa9rz7Q9naz9WLc3jnUSdx1Ba+p6I1iPl0fyBf/2cR1tHQmzETH7kn48Etayvm1yX2GXfsNeZERUBBwxIyIiIioIXpgRERERFQSjTDoXihZ99vZTA/oLDMLbDCo06Nse92Tr7o9/LvXvA4CSRpT2Jf5Io0kXLCOXXM/DJj7vaEGArduyHU/3pOstnQbK4sa++FGnkWolY8fwON/pdCNK29bOF1EysiSiIuKIGREREVFB8MKMiIiIqCAYZdKFMe6IKu3+0vqpAdnxZ9/tbD0jCs08PirFcWh4TN5q06z9eWTGgu2TY8TuspO63YQ9xnrv57SRJCNLIioijpgRERERFQRHzIgyPIgRlUH3OWjEre+4k+4ja39puBGxcLTOdIY4P+FIV7z9AY52cUSMiB5FHDEjIiIiKghemBEREREVBKNMojF6WPHZMI+T1Y8tr3G+ogc1iT0R0XnBETMiIiKiguCIGVGBPMiCg9O0wRj3cyAiopNxxIyIiIioIHhhRkRERFQQjDKJLgjGiURExccRMyIiIqKC4IUZERERUUHwwoyIiIioIHhhRkRERFQQvDAjIiIiKghemBEREREVBC/MiIiIiAqCF2ZEREREBcELMyIiIqKC4IUZERERUUHwwoyIiIioIHhhRkRERFQQvDAjIiIiKghemBEREREVBC/MiIiIiAqCF2ZEREREBcELMyIiIqKC4IUZERERUUHwwoyIiIioIHhhRkRERFQQvDAjIiIiKghemBEREREVBC/MiIiIiAqCF2ZEREREBTHShZlzbsk599vOudd0uXjCsSXn3B84535jlMckIiIiOq9GHTH7EoDf8d4/C+B3dD3LLwB4ZcTHIyIiIjq3Rr0wexHAr+jPvwLgr6Yd5Jy7BuC/BPDVER+PiIiI6Nwa9cLskvf+FgDoci3juH8I4H8C0Bnx8YiIiIjOrfKgA5xz/xbA5ZRdfz/PAzjn/isAG9777zrnfjrH8S8BeAkAUJnO8xBERERE58LACzPv/V/K2uecu+OcW/fe33LOrQPYSDnspwD81865zwGYADDnnPvn3vv/JuPxXgbwMgBEUys+z4sgIiIiOg9GjTK/AeDn9OefA/Dr4QHe+//Ze3/Ne/8EgM8D+HdZF2VEREREF9moF2ZfBvBp59xrAD6t63DOXXHOfXPUJ0dERER0kQyMMk/ivd8E8BdTtr8P4HMp2/8DgP8wymMSERERnVfs/E9ERERUELwwIyIiIioIXpgRERERFQQvzIiIiIgKghdmRERERAXBCzMiIiKiguCFGREREVFB8MKMiIiIqCB4YUZERERUELwwIyIiIioIXpgRERERFQQvzIiIiIgKghdmRERERAXBCzMiIiKiguCFGREREVFB8MKMiIiIqCB4YUZERERUELwwIyIiIioIXpgRERERFQQvzIiIiIgKghdmRERERAXBCzMiIiKignDe+7N+Dpmcc3sAfnzWz6NgVgDcO+snUUA8L+l4XtLxvPTjOUnH85KO5yXdB7z3s6PcQXlcz+QB+bH3/pNn/SSKxDn3HZ6Tfjwv6Xhe0vG89OM5Scfzko7nJZ1z7juj3gejTCIiIqKC4IUZERERUUEU/cLs5bN+AgXEc5KO5yUdz0s6npd+PCfpeF7S8bykG/m8FPrL/0REREQXSdFHzIiIiIgujDO5MHPOfcY592Pn3OvOuS+l7HfOuf9D9/+Rc+7jeW/7KMtxXv6mno8/cs79rnPuJ3r2ve2c+4Fz7g/HURVSJDnOy08753b0tf+hc+6X8t72UZXjnPyPPefjj51zbefcku47z++VrznnNpxzf5yx/8J9tuQ4Jxf1c2XQeblwnytArvNy4T5bnHOPOef+vXPuFefcD51zv5ByzPg+W7z3D/U/ACUAbwB4CkAVwPcBvBAc8zkAvwnAAfgzAH4v720f1f9ynpf/DMCi/vxZOy+6/jaAlbN+HWd0Xn4awG+c5raP4n/Dvi4AfwXAvzvv7xV9bf85gI8D+OOM/Rfxs2XQOblwnys5z8uF+lzJe16CYy/EZwuAdQAf159nAbz6IK9bzmLE7FMAXvfev+m9bwD4OoAXg2NeBPDPvPgWgAXn3HrO2z6qBr427/3veu+3dPVbAK495Od4Fkb5nZ/X98uwr+tnAPzLh/LMzpj3/v8BcP+EQy7cZ8ugc3JBP1fyvFeynNv3CjD0ebkQny3e+1ve++/pz3sAXgFwNThsbJ8tZ3FhdhXAuz3rN9H/ArOOyXPbR9Wwr+3vQK7OjQfwb5xz33XOvfQAnt9ZyXte/qxz7vvOud90zn1wyNs+anK/LufcFIDPAPjVns3n9b2Sx0X8bBnGRflcyesifa4M5aJ+tjjnngDwMQC/F+wa22fLWXT+dynbwtLQrGPy3PZRlfu1Oef+AuQD9M/1bP4p7/37zrk1AL/tnPsT/ZfPoy7PefkegMe99/vOuc8B+NcAns1520fRMK/rrwD4T9773n8Bn9f3Sh4X8bMllwv2uZLHRftcGdaF+2xxzs1ALkT/nvd+N9ydcpNTfbacxYjZTQCP9axfA/B+zmPy3PZRleu1Oec+AuCrAF703m/adu/9+7rcAPCvIMOn58HA8+K93/Xe7+vP3wRQcc6t5LntI2qY1/V5BFHDOX6v5HERP1sGuoCfKwNdwM+VYV2ozxbnXAVyUfYvvPe/lnLI2D5bzuLC7NsAnnXOPemcq0J+ud8IjvkGgP9Wqxz+DIAd7/2tnLd9VA18bc656wB+DcAXvPev9myfds7N2s8A/jKA1IqaR1Ce83LZOef0509B3tebeW77iMr1upxz8wD+PIBf79l2nt8reVzEz5YTXdDPlYEu4OdKbhfts0XfB/8EwCve+3+QcdjYPlseepTpvW85574I4Lcg1Qpf897/0Dn387r/KwC+CalweB3AIYC/fdJtH/ZreBBynpdfArAM4B/r50XLyySylwD8K91WBvB/ee//7zN4GWOX87z8DQB/1znXAlAH8Hkv5TDn8v2S85wAwF8D8G+89wc9Nz+37xUAcM79S0g13Ypz7iaA/xVABbi4ny05zsmF+1wBcp2XC/W5YnKcF+Difbb8FIAvAPiBc+4Pddv/AuA6MP7PFnb+JyIiIioIdv4nIiIiKghemBEREREVBC/MiIiIiAqCF2ZEREREBcELMyIiIqKC4IUZERERUUHwwoyIiIioIHhhRkRERFQQ/z/4jQgaVzxmjQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make plot of scalar field\n",
    "x = domain.grid(0,scales=domain.dealias)\n",
    "y = domain.grid(1,scales=domain.dealias)\n",
    "xm, ym = np.meshgrid(x,y)\n",
    "fig, axis = plt.subplots(figsize=(10,5))\n",
    "p = axis.pcolormesh(xm, ym, S['g'].T, cmap='RdBu_r');\n",
    "axis.set_xlim([0,2.])\n",
    "axis.set_ylim([-0.5,0.5])\n",
    "\n",
    "logger.info('Starting loop')\n",
    "start_time = time.time()\n",
    "while solver.ok:\n",
    "    dt = cfl.compute_dt()\n",
    "    solver.step(dt)\n",
    "    if solver.iteration % 10 == 0:\n",
    "        # Update plot of scalar field\n",
    "        p.set_array(np.ravel(S['g'][:-1,:-1].T))\n",
    "        display.clear_output()\n",
    "        display.display(plt.gcf())\n",
    "        logger.info('Iteration: %i, Time: %e, dt: %e' %(solver.iteration, solver.sim_time, dt))\n",
    "\n",
    "end_time = time.time()\n",
    "\n",
    "p.set_array(np.ravel(S['g'][:-1,:-1].T))\n",
    "display.clear_output()\n",
    "# Print statistics\n",
    "logger.info('Run time: %f' %(end_time-start_time))\n",
    "logger.info('Iterations: %i' %solver.iteration)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "As an example of doing some analysis, we will load in the horizontally averaged profiles of the scalar field $s$ and plot them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in the data\n",
    "f = h5py.File('analysis_tasks/analysis_tasks_s1/analysis_tasks_s1_p0.h5','r')\n",
    "y = f['/scales/y/1.0'][:]\n",
    "t = f['scales']['sim_time'][:]\n",
    "S_ave = f['tasks']['S profile'][:]\n",
    "f.close()\n",
    "\n",
    "S_ave = S_ave[:,0,:] # remove length-one x dimension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEvCAYAAAAD0BVUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABO40lEQVR4nO3deXzcV33v/9eZGWlmtIz2fZcsW7It2/Ee23HsJCaJk5CFJU4CCZA0BEpLufdXSEOhpJcSUm5boMCFEOhN4QZKIU0CDUnI6sRrvMu7tVubtS+jGc16fn98xyPJq/aR7M/z8dBjtvP9zpmvrXnrfM/5nqO01gghhBDTzRTpCgghhLg6SQAJIYSICAkgIYQQESEBJIQQIiIkgIQQQkSEBJAQQoiIsES6AjNJamqqLiwsjHQ1hBBiVug83URdW2uH1jptPNtLAA1TWFjInj17Il0NIYSYFZ770hN86rtP1Y93ezkFJ4QQIiIkgIQQQkSEBJAQQoiIkAASQggRERJAQgghIkICSAghRERIAAkhhIgICSAhhBARIQEkhBAiIiSAhBBCRIQEkBBCiIiQABJCCBEREkBCCCEiQgJICCFEREgACSGEiAgJICGEEBExKwNIKXWLUuqEUqpKKfX4JcqtUEoFlFIfnc76CSGEuLxZF0BKKTPwQ+BWYD5wn1Jq/kXKPQ28Nr01FEIIMRqzLoCAlUCV1rpGa+0Ffg3ceYFyfwH8DmibzsoJIYQYndkYQDnA6WGPG0PPhSmlcoC7gR9fbmdKqUeVUnuUUnva29sntaJCCCEubjYGkLrAc/qcx98FvqK1DlxuZ1rrZ7TWy7XWy9PS0iajfkIIIUbBEukKjEMjkDfscS7QfE6Z5cCvlVIAqcBmpZRfa/3itNRQCCHEZc3GAPoAKFVKFQFNwBbg/uEFtNZFZ+8rpf4v8AcJHyGEmFlmXQBprf1KqS9gjG4zAz/XWh9RSj0Wev2y/T5CCCEib9YFEIDW+hXglXOeu2DwaK0/NR11EkIIMTazcRCCEEKIK4AEkBBCiIiQABJCCBEREkBCCCHG6UKXZY6eBJAQQoiIkAASQggRERJAQgghIkICSAghRERIAAkhhIgICSAhhBARIQEkhBBinGQYthBCiIiQABJCCBEB564EOlYSQEIIIcZFaWkBCSGEiAA9wQiRABJCCDE+WgJICCFEBOgJrmkqASSEEGJcdDB6QttLAAkhhBiXoLZNaHsJICGEEOMS1HET2l4CSAghxJgFA0H8gcwJ7UMCSAghxJh1NDrRSB+QEEKIaRQMBNnxX9VAYEL7kQASQggxJu//ZxWNx7txWN6Y0H4kgIQQQoxKIBBk/+sNVL7TyOKb8ohRlRPa38SuIhJCCHHF6213c3RbM8e3t+Dq81JYkcKae+bw4lvBCe1XAkgIIcR5Av4gNQfaOfp+M43Hu1EKCipSmb8um4IFyZhMCpNPAkgIIcQk8HsDdDQ6qdnfzvGdLbj7fcQn21h5RxHla7KISxp54anJP7FBCBJAQghxFTobNu0N/bTV99He0E9Xiwsd1JhMisLFRmsnr9xo7VyIxSsBJIQQ4hIuFTYA9vgo0vIdFC1OIy0/nsziBGIcl7/GJ2rAO6F6SQAJIcQVxDvop6t54DJhEx8Om7T8eOKSrCg19sXlrAO+CdVVAkgIIWahgC9I9xkXXc1OOpsH6GoeoKvZSV/HYLjM2bApXJRKeoFjQmFzIfZ+aQEJIcQVKxjU9LW76WoeoLPZGbodoPeMi2CoVWMyKRIzY0gvdFC+Jovk7LhJD5tz+bu7ieuRABJCiFlPa42z2zMiaLqaB+hqGSBwdrizAkeqnZTsWIqXpJKSHUdydiyJGTGYLdM7r8Dg4cMT3sesDCCl1C3A9wAz8KzW+tvnvP4A8JXQQyfwOa31wemtpRBCXNjZsGmr76Otvp/20K3H5Q+XiU2IJiUnjorrc0jOjiMlJ5akzFiirOYI1nyIa9cu9AT3MesCSCllBn4IbAIagQ+UUi9rrY8OK1YLXK+17lZK3Qo8A6ya/toKIQQM9AyFTVt9P+0Nfbj7jQ58k0mRkhtHybJ0UnPiSMkxWjW22KgI1/ri/F1ddD//K6rmxcCJ8e9n1gUQsBKo0lrXACilfg3cCYQDSGu9fVj5nUDutNZQCHHVcvV5w6PPjMDpw9Vr9JUokyI5K5bCilTSC+JJK3CQkhOLJWpmtGpGq+PHPybo8fDm+mR4efz7mY0BlAOcHva4kUu3bh4G/jilNRJCXJV0UNPVMkDTyR6aT3VzprYPZ7fHeFFBUkYMeWXJpBXEk17gIDUvjqjo2RU253IfPEjPr35N4j1305j89oT2NRsD6EJDOi54KlIptREjgNZddGdKPQo8CpCfnz8Z9RNCXKGCQU1no5Omk900n+qhuaoHz4DRbxOXbCWrJIH0QgfpBfGk5sUTbZuNX7EXprWm59e/pvVbTxGVnk7qF/6Cnlcm9rf9bDw6jUDesMe5QPO5hZRSi4BngVu11p0X25nW+hmMPiKWL18+0T41IcQVJBgI0t7gpOmUETgtVb143UbgOFJtFC1OI6c0kezSRByp9gjXduoE3W5a/u7v6Hv598Rev56cp58GRzyd6uobhv0BUKqUKgKagC3A/cMLKKXygReAT2qtT05/FYUQs5EOajoanTQc7aT5ZA8t1b34PMZ8Z4kZMcxZnk72nERy5iaeNzHnlUgHg7j27OHMN/8Bz6lTpP7lX5D62GMok4nWgVYCamJ/s8+6ANJa+5VSXwBewxiG/XOt9RGl1GOh138MfB1IAX4UugjLr7VeHqk6CyFmrsEBH6ePdtFwpJP6o124+4y/6pOzY5m3OpPsUAsnNsEa4ZpOH09tLb0vv0zfSy/ja27GnJRE3k9/Sty6teEyJ7sn/rf9rAsgAK31K8Ar5zz342H3HwEeme56CSFmPh3UtNX303C0k/rDnbTV9aE1WGMt5Jcnk78whbzy5KsqcMCY2aDvlVfofellBg8dApOJ2DVrSPvSXxF/442YYmJGlN/VsouoCXZazMoAEkKIsXD3e2kItXIajnYx6PSBgvQCB8s3F5K/IIX0QsdFlx24UgXdbpxb36P3pZdwbt0Kfj/WsjLSv/xlHLffRlR6+kW33dWyi4WeAPsn8P4SQEKIK1Jfh5uaA+3UHGintboXrY3JOfMXJFOwIIW8+cnY4y6/5MCVQmuNr7ER94GDuA8cwH3gAIMnToDfjyUtjeQHHyThzg9jmzfvsvtq6GvgRPcJ/tw1wC8mUCcJICHEFUFrTWfTQDh0OhudAKTkxLJscyFFi1JJy4tHXSWtnKDbzeDhw7gOHDBC5+BBAh0dAKiYGOwVFaQ8/DCxq1YSs2oVyjz665N+sP8H2M1WPurs43MTqKMEkBBi1goGNa01vdQcaKf2QLuxFIGCrOIE1nxkDsVLUklIi7n8jmY5rTW+pibc+w+c17oBiCrIJ27tGuxLlmBfsgRraSnKMr6v/6OdR/lj3R/5s5wbSK06NaF6SwAJIWYVHdS0VPdw8oM2ava34e73YTIrcsuSWXpzAYWLUq/oAQQ6EMDb0IDnxAkGjx/Hc/wE7sOHh1o3drvRuvnMZ0KBsxhLcvKkvLcn4OHp3U+TaE3kPpU14f1JAAkhZjytNR2nnZzc3UrV3jac3R4sUSYKKlIpuSaNgoUpRNuvvK+zgNM5ImgGT5zAc+oU2u02CpjNRBcVTlrr5lJcPhd/+dZfsq9tH99c+02i3/4ZtcEMoG/c+7zy/sWEEFeM7tYBTn1whlN72ug548JkUuQvSObau0soXJR6xUx1o4NBfI2NRtCcOMngCSNwfI2N4TKmhARs8+aR+LGPYps3D+u8MqylczBZp7611+vp5fNvfJ4jnUf4h3X/wIdzNxJofpjfBDcB4z8Nd2X86wkhrhgDPR5O7j7DyQ9a6TjtBAU5cxNZclMeJUvTZ/QyBaMRdLvxnDzJ4LHj4aDxnDhB0OUyCihFdGEhtoULSfzoR7DOm4etrAxLZuaUrW56Kc3OZr7w1heo663jnzb8Ezfm3whHXsSsfWy3rAJ+Pe59SwAJISIu4A9SX9nJse3N1B/uRGtIL3Sw7mOlzFmWTmzi7OzT8Xd0GEFz/BieY8cZPH4cb10dBI0VTk1xcVjnzSPhrruwlhlBYy0txWSP/LxyLc4Wnq18lheqXiDKFMUPb/wh12ZfC1qjP/gpPSqBQPbKCb2HBJAQImI6m50c297CyV2tuPt9xCZEc83NBZRfm0VixuwZvaYDAbz1DXiOHwsFjhE6gfaOcJmo7Gys5eU4br0VW3kZ1rIyonJyItKquZQmZxPPVj7Li1UvAnDPnHt4pOIRsuJCgw6O/R5V9z7/5Ps0H6rIkeuAhBCzh8ftp2rPGY5tb+FMbR8ms6JoUSpla7LIn5+MyWyKdBUvKehyGafQjh8fat2cHDYwICoKa0kJcWvXhYKmHFvZPMwJCZGt+GWc7j/Nzyp/xktVL6GU4iOlH+GRikfIjM0cKuQbhNe/Spu9hP/03siOiomNhJMAEkJMi45GJ5Vvn+bk7jP4fUGSs2NZ+9E5zFuViT1+Zs5IEPR68Rw/jruyksHKw7grK/HW1IA2JkEzORzYyspI+vjHjKApL8NaXIyKnpmf51wNfQ1sbdzK1satfND6ASZl4mPzPsZnFn5mZPCA8Zlf/1voaeDvLd9gTWkGybET+5wSQEKIKRMMBKk92MGhtxtpPtWDJcrE3JUZzF+XQ3ph/Iw6/aQDAbw1NbgrD+OuPMRg5WHjYk6fDwBzSgr2igoct9yCbX65MTAgO3tGfYbL8QV87Dmzh/ea3uO9xveo66sDoCihiE8u+CQPlD1ARmzG+RsG/PCHL8L+X7IrYwt/qJ/L8+uLJ1wfCSAhxKRzO70cfb+Zw+824ez2EJ9iY809cyhfmzUjRrEZMwc0M1h5CHflYQYrKxk8ciQ8Es0UG4tt4UJSHnoQW8Ui7BULsWRlzaqwOavd1c57Te+xtXErO5p34PK7iDZFsyJzBVvKtrA+dz158XkX34FvEH73MBz/Ay3XfIn7di7nvpX5rClJnXDdJICEEJOms8nJgTdPc2r3GQL+ILllSVx371wKF6VGdKbpgNNpzIe2f3+4dRPo7gZARUVhLS8n4e67sVUsxF5RQXRREco0s/uiLqbT3cmRziMcbD/Ie43vcazrGAAZMRncVnwb63PXszJzJTFRoxjk0d8KLzwKte/i3fQUD+1eQHq8n7/ZXDYpdZUAEkJMWEt1L/teraOushNLtImyNVlUbMghJTsuIvXxtbbi2rsX9779uPbtw3PihDH02WQyBghs3Ii9YiG2ikXY5pbOmj6bc3UNdnG08yhHO49ypOMIRzqPcMZ1BgCTMrEkbQlfXPpFrsu5jrlJc0ffghvogPf/BT54FoJ+em7+Vz6xp4hTbX38/KEVOGyT04qVABJCjIvWmtNHu9j7aj3Np3qwxUax8o4iKjbkTutpNh0I4Kmqwr1vH669+3Dt24u/uQUIzfq8eBGpn/sc9qXXYF+8BHNc7LTVbTL1DPYYQdN5JHzbMtASfr3QUcjSjKUsSFnA/JT5lCeXExc9xj8AXF2w/V9h10/A74ZF93K49DE+9WIHgz4Xzz64nI1lF18jaKwkgIQQYxIMamr2t7P31To6TjuJS7Ky7mOlzF+XTZR19FP6j/v93W7chypx7zcCx33gAMH+fgAsaWnYly0j5lOfwr50GbayeVMyL9pU6/X0hoPm7E+Tsyn8en58PkvSlnB/2f0sSF1AWXIZ8dHx43/DwV7Y8SPY+SPw9MPCe+D6x/ldQwx/86tKMhNsPP9nq5ibMYH3uIDZ9y8jhIgIHdRU7Wtj9+9r6TnjIjEjho2fLGPeqkzMlqnrLwl6PLj3H2Bg105cu3bjrqwMj0yLnlOC49ZbiVm2FPvSpUTl5s6qgQK9nl5qe2up7qmmureamp4aqnuraR1oDZfJjctlYepCPj7v4yxIWUB5SjmOaMfE33ygA079CU6+ClVvgrcfym5Hb/gbdruy+Okrtbxx7BTXFqfwoweWkjTBIdcXIgEkhLgkrTV1hzrY9XItnU1OkrNjufnPFlJ8TdqUDCzQXi/uykoGdu0yAmf/frTXCyYTtgULSHnoQaOVc801mBMTJ/39p0LXYBc1PTXU9NaMCJt2d3u4jM1soyihiOUZy5mTOIf5KfOZnzKfBOskXcCqNbQdhRN/hJOvQeMHgIa4DFhwF/5lD/NKZwbP/raGQ431JMVE8aWb5vL5jSVETdHFwRJAQogL0lpz+lgXu16qoa2+n4R0O5sens+cZRmTGjza72fw6FEjcHbuwrVvnzGrgFJYy8pIuu8+YlavImb5cszxk3sKaDJprelwd1DdW011z1Brpqanhm5Pd7hcjCWGksQS1mSvoSSxhJLEEooTismOy8akJvmL3jcIde8ZrZyTr0HvaeP57Gtgw+Mw92b6k+bzH3ua+Ldf1tHU00xxaizfvGshH1maiz16ak+pSgAJIc7TXNXDzheraanqJS7ZysZPllG2OnNSpsnRWuOtqsK5bZsROHv2EHQay2dbS+eQeM89xKxaScyKFViSkib8fpPNH/TT7Gymrq+O2t7aEafQ+r394XLx0fHMSZzDDfk3UJxQHA6bjJiMqTtN6HHCmSPQegiq34Kad8DngqgYKLkBrv8ylH4Ijz2NAw09/GnfGf7jg3fo9/hZWZTMNz68gBvL0qdtyLwEkBAirKtlgB3/VU3doQ5iEqJZv2Uu89dmY46aWPD4u7tx7diB8/1tDGzbhv+MMVQ4uqAAx223EbtqJTErV2JJnfjFjZOlZ7AnHDJ1fXXU9dZR11dHQ38D/qA/XC7ZlkxxQjGbizaPCJoUW8rUBY3W0N8CrZVG2LRWQuth6KoBjGmCSMiHaz4Bc2/Gm7uGg62D7KjuZOev69hbvx+PP4jZpNhckcWfXVfEotzEqanrJUgACSEY6PGw+/c1HNveQpTVzKo7i1l8Yx5R4zwFo30+3AcO4Ny2jYH3tzF45AhojSkhgdhrryV27Rri1q4lKjt7kj/J2PgCPhr6G8LhMjxoejw94XIWk4X8+HwKHYVsyNtAoaOQooQiChwFJNmmuJUW8EHHSSNgwmFTCe6uoTJJRZBZAYvvg8wKvKnzOdQXy87aLna+08We+ncZ9AVRCsozHTywqoBrS1JYWZhMQkzkZqaQABLiKuZ1+9n3ej0H3zhNMKip2JjL8s2F2OPGPuLJW18fDhzXrl0EBwbAbMa+eDGpf/EF4tauxbZwIco89UO1hzvbN3Oh1kyTs4mgDobLptpTKXQUclPBTeGQKXQUkh2XjcU0DV+Xg72hoKmEM6GgaTsGAa/xusUG6fOh/HbIXGSETvp8fFFxHGrsZWdNJzvf72RP3VHcvgAAZZnxbFmRz7UlKawqSiYxZuZcdCsBJMRVKBAIcmRrMx/8dy2DTh+lKzJY9eFiEtJGvxBa0OPBtfsDnO++i3PrVnwNDQBE5ebiuON2YteuJXb16mkbOOD2u6nvqx8RMHW9ddT31eP0OcPlbGYbBY4C5qfMZ3PRZgoTCilyGK2ZMV+4OVbBIDhbobseuuugp37k/b6ha32ISYWsRbDqsXDY9MUVUNvpobZjgJqOAWqrB6jtOEh120A4cOZlxHPvijxWFyezsihlwjNWTyUJICGuIlprag92sP2FKnrb3OTMTWTNR+aQXjC660p8zc04t27F+e5WBnbuRLvdKJuN2FWrSH7oQeLWrSMqP3/K+j6COkjLQAv1vfXU9tWOOHU2/NoZgKzYLAodhdxRcgeFjsJw0GTEZkz+aLPh3D2hYKk7P2h6GiDgGVZYgSMbEgug6HpInYMvbQGN1jmcHIilttNFbfsAtTUD1HQ00uGsDm9pUpCbFENRaiwrViazojCZVUXJpMTNntVjJYCEuEqcqe1j2+9O0VLVS1JmDLf9+SIKFl66o1z7/UZfzrvv4nx3K56TJwGjlZN4zz3EXb+emJUrMdlsk1rXPm8f9b31I0+b9dXR0NeAZ9gXeFxUHIWOQpZnLA+HTKGjkAJHATbL5NYpzO+BntOhYKkbCpqzoTPYO7K8LRGSCiC9HObdCkkFBBMKaY/KpMqTRHW3j5r2AWo7Bqg9NUBjt4ugPhHePDXOSnFqLDeWpVOUFktRaiwlabHkJcdgtUzv6czJpnRoYSUBy5cv13v27Il0NYSYVH0dbna+WM2pPW3Y46NYeUcx89dmXXRItb+7m4GtW43QeX8bwb4+sFiIWbaMuOuvJ+769UQXF0+4leMP+mlyNlHXOxQyZ2+7Boc62M3KTG58rhEwjkIKEgrC/TOTPtJMaxjsMWaB7m+Bvhaj1TL8VFl/C+GRZgDmaKMFk1QASYWQWEAwMZ9eay6t5nRavTba+zzUdYZCJvTj8Q/1PcVEmylKNcKlOC2O4tD9wtRYEuyRX77iUpRSe7XWy8ezrbSAhLhCDQ742PPHOirfacSkFMs3F3LNh/KJtp3/a++prcX51tv0v/0W7n37IRjEnJpK/E03Ebd+PbFr14y7L8fpdY64ZubsT31//YjhzEnWJAoTCrk+9/pwS6YwoZC8uDyizJPwJex1GeHR3zIUMP2t0Nc88rHffc6GQ6fJAoXrccbk0h2dRas5k0bSafDE0+b00dbvoa1mkPZ+Dx1OL4HgGeBMeC8WkyI/2Thltm5OKkVpsRSnxlGcFkt6vHVWTSE0WSSAhLjC+H0BKt9uYu+rdXjcfsqvzWLlHcXEJQ31DZw9tdb/9ts433obb20tANayMlIf+yxxGzdiW7Bg1GviBHWQ1oFWozXTNzJohk83Y1Zm8uLzKEooYn3eeoocRRQlGD/jnnIm4AsFSOuwgGkZ2YrpbwVP7/nbWuxoRxaB2AzcqYvpz7qRLlMybTqJ5kAi9T4HJwcdNPdr2hs99Jz0nbODdpRqJyXWSnq8lbR4K+WZDtIdVtLjbeHn0uNtZCXapmxKm9lKAkiIK4QOak7tOcPOF2vo7xokf0EKa+4pISXHGNkVcA4w8P77ON9+G+e77xLo6YGoKGJXriTpgQeI37iBqJycS77H2ZFmtb214VNnZwcDDAYGw+Xio+MpSihiTfaacMCMuTUTDIKrY1iItFw4ZAbaz9/WZEHHZ+KPycAdV0xf0kq6zCnhYKnzxlPldlA/YKat3Yu7OXDeLqLNJiM8HGaKUq2sKk4Oh0q6w0panI10h5WU2GgsEizjMqoAUkot1lofnOrKCCHGp/F4F9tfqKa9oZ/UvDg2PriEvLJkfG1tdP/6v+l/801cO3eifT7MCQnEXr+e+BtuIHbdOsxx5w897hrsora3lpreGmp6asKtmeaB5nAZhSI7LpuihCJWZK4wgibUokm2JV/4lJLWxpozzjYYaAvdtg973D70vPMMDDtFB6BR+O0puK3pOK1p9CbMpSspmTaSaAkm0ehPoNbjoM5lo73dhz94fh93vNVCmsNosSzKCwVKKFTS422hFouVBHvUVXlabDqNahCCUsoH/BvwNa31mcuVn61kEIKYbTqbnGx/oZqGI53EJVlZfWcxBWkunG++Rf+bbzB48BAAUfn5xN9wA/E33oD9mmtQFkt4SPPZWZrPhkxNb82IWQDsFnu403/4T358vjHSLBgEV+flA2WgHT3QjjonVACCyozLkkS/JYkeUyIdJNIaTKLRn0idz0Gdx8EZnUQHCfjP+bvZpCAxJprEmCgS7VEkxUSTFBs9LFhGngqb6gk2rzYTGYQw2gAaAGyAE3gK+BettefSW80+EkBitujvGmT372s4vrMVq83CosVR5LXvwPXWG+H+HNvChcTfdCO2jetpTrOE+2bOhs25p82SrEnhcCl2FFJsTaHYHENmIIBpYGTAaGcbQafx2OTuRA2bTeAsPxb6zEl0qQQ6dAJngg6affG0a+NxO8Zth3bQQxwaEw6bhcSYaJJiosKhkhQTTYI9iqSYKJJiz94PBU5MNPFWy7RNninONx0BlA38I3Bf6Kl64HGt9W/G86YTpZS6BfgeYAae1Vp/+5zXVej1zYAL+JTWet/l9isBJGY6j8vHvtfqOfhWIzoQpMR2mpx9/w9T62mwmAles4D25YUcLY/lmLmd2r5aGvsbCeihPo5sWypFtlSKLfEUEUW+L0iey43D2YlpoAOLu50obw+K878bPETTSQLtQQdtoSDpCAfJ2ccOBqKSMduTSIwdCorEc4LDCJmz96Nx2CzSlzILTXkADXujlcB3gdUYA+F3Al/SWu8ez5uPh1LKDJwENgGNwAfAfVrro8PKbAb+AiOAVgHf01qvuty+JYDETBXwBTn4ejV7X23A64OMjr0UV72ERfdRUxbHthI/WwvcuGxGS8CCIjsYRb5PU+z1UjroZJ7XSYHPT8w5v/MD2np+kJBAu06gx5SI15qCz5aGjk3FGptAUqx1WOskigR7dLh1kmiPIiEmatZfIClGb9quAwoFzRql1P3At4FrgR1KqeeBv9FaN46nEmO0EqjSWtcAKKV+DdwJHB1W5k7g37WRrjuVUolKqSytdcs01E+IMfP4A/S5/fT19eLqbmWwp5nWniramo9CZTR+9zr8lmSSu46R3vwSR3ObePkORXW+Jk93Uuz182cuH8W9Pgq8fmJ9dnp0Ap3aQZfKoN+cxLaoJN6KS8ZnTSEYkwqxaZjjM4iJc+CwW0iwR+GwR5E7rKVijzJLR7yYMuMahq21fl4p9QLwFeCvgQeAe5RS/xt4WmvtmsQ6nisHOD3scSNGK+dyZXKA8wJIKfUo8ChAfn7+pFZUXD201rh9AXrdPnrdPvrcfuN2wI2nr51Afxva2Y7J1Y7Z3YnV24Xd24XV34Xb3ENflJOOaB+N0YpeZxTptWbKm8owJ9zFQFwuNm8DyvtL6ktrqVkTR4qphNvMmcSbstH2NFRcOhZHOmZHBp0JqfhjbSTaoyiwR2GLMkmIiBlp3NcBaa0HgSeVUj8FngbuB/4WeEQp9YTW+rlJquO5LvSbdO55xNGUMZ7U+hngGTBOwU2samK201rT7/HTM+Cj2+Wlx+2jx+WlxxV67PLhdPYTdLZjdrUTPdiBzdNJrL+LJN1LmuolhT5SVC8lqo8knJiUZkAp6qKiqIm2UBMVRU1sNFVJVposJgJoilujWVFpYdMphcOdQ1XxnXTllGM1D7B8vYUVH30Qk/lTkT48QkyqCV+IqrVuBj6plPoh8CugAPi5UuovgP+htd460fc4RyOQN+xxLtA8jjLiCna2RdLtGgqQoRA5Gyg+et1eul0+3AN9WFwdRHs6SNY9pKo+UuklVRk/5aqXVHpJM/USz7lTtQBmGLTE0WpLoTrGwVFbDqej82kw+ajXLtr8Q8sBWJSZgtg8NpxJZNHBQbL2nyaqoxd3TBr1Sz/JCUsJVruZdbcXs3B9zoRXIxViphp3ACmlSoDlwIrQzzVA7NmXgaXA20qp/8QIoskKgA+AUqVUEdAEbMFofQ33MvCFUP/QKqBX+n9mL48/QG8oMHpc3qFQcRuB0hsKlm6XL3y/x+0l2j9gBAi9RqAoo4WSSi9l5j4yTMZzSboHmw4NRz7nIn2/NZFgjHGKyxw/H1N8BoGYFJqtNmoJUBMcoMbTTa2rlZq+Ovq8fUAv0Itd2SlyFLEyodhYqjkqm7yjHdi2H8K19T2C/VUom43odTdQm3sTJ0/bUCbF0hvzWHpzAVa7TFQirmyjnQkhh6GgWR76SRxeJHTrBvZhjI5zA38OfBz4kFLqYa31f020wlprv1LqC8BrGMOwf661PqKUeiz0+o+BVzBGwFVhDMP+9ETfV0wOp8dPe7+HrgEP3QO+C57i6nF76R4w+lK6XV5c3qEhxLG4yVDdpIVaJ5nmPuZG9ZNl7ifN1EsKvSTSQ1xUN1GW8y9V0yiISUHFpUNsJsQtgth0iEsL3aZDbBrEpTMYHUe9q2XEBZo1vVXUn35zxJIAybZkihOKubnwZopDYVOcWEx6TDrBjk76336b/l++iWv7Dvw+H+6kJOI3bcK+4UZqvPm892Yzvno/ZddmsvKOIuKSpmgZASFmmNFeBzR8oqTh/SvVGGFz9ueg1to/bLsE4FvA54AgcLvW+tVJqPeUkGHY4xMMarpdXmM24H4P7f0e2voHaesbdj/0/PAwGc5u8lNic1Js7SM/uo8cc48RNLqLpGAn8b4OYjztRPkHzt9YmYzQGBEk5wcKsekQkwLmob+7tNZ0DnaOWA4gPO2Msxkd6jo8O+3M2YApSiiiONG4P3wSTa01npOncL79Fv1vvz00E0FeHvE33mhcGLp4CSf3trPrpRqc3R4KKlK49q6hOduEmE2m40LUINCHcforHDha685RVvAvMa4f2qa1vm48FZ0OEkAjef1BOpxGqLT1DdLu9NDWdzZkhkKlvd9z0Tm3MuItzIlxU2zvp8DSS7a5hzS6SPB3Eudtx+puw+I6g8nddX4FzNEQnwnx2aHbrGG3GUMBY0+Gy8za7Av6ON1/emgCzdAkmrW9tfR7+8PlbGZbeCmAooQiCh2FFCcWU+AowG658HLV2uvFtXcv/W+9jfOtt/A1Gcsq2xYtIn7jBuJuuAHr3LkopTh9rIvtL1TRcdpJWn48az4yh9x5SaP/RxFihpmOAFoAHNUTWL1OKdULoLUe55zrU+9qCaABjz8cKkMtFqOl0t5vhEy700PXgPeC26fERpMWF01xnJcSWz/5Ub1kmXrCrZU4bzu2wXZM/S3G9C3nTtOiTBCXMSxMskaGiyP02J4EYxw+3OvpvWBrprG/Ef9Q45w0e1o4YIbPb5YZmzmq5ZoDvb04t76H8+23cG59j6DTibJaiV2zhriNG4jbsIGo9PRweWPOtioajnQRn2Jj9V3FlC7LQMkUMmKWm/ILUbXWR8az83N0M3JkmpgCWmu6BrzUdgxQE1p5saHLRXvfUMAMXOA0WJRZkRZnJc1hIz8lhtX5MRRHdVGg2sgMniHF30L8YCvR7jOYnKG1V3ovEFD2ZGPxrvhMyFxw4YCJSwfT+K+UH76S5vBVNGt7a0espGkxWSiIL2BO4hxuKrgpPFtzYUIh8dFjW1xNa423thbnu1txvvMOrj17IBAwFm275WZjZulrr8VkH9lKcnZ7jDnbdrQQbbew5iNzWLQhV0a2CcH0rgf0KMYUPmISuLz+oeV9242wqekYoLbdSd/g0F/6FpMiJ8lORryNBTkJoRmCbaTHWci19JAVPEOyr4XYgUbU2WWH2+rA2TryDS02SMg1AiRv9VArZXgrJi4DoianA11rTbenOxwy4du+Ok73nx6xkmaiNZGihCI25G0IB0xRQhE5cTlYTOP/Lx70eHDt3o3znXdxbt2K77RxbbO1tJSURx4h/oaN2CoqLrhom3fQz/7XGzjwpwaCWrPoxjyW31qILXZmL68sxHQa01xwV7qZdgrOFwjS2O2mtsNJTfvQWvI17QO09g2OKJuVYKM4LTa0rnwcxSkxlMT7yKYNS++w9ezPhkxPAwSHr+6ojIA5Z217kgqNx7Hpl+1nGQ9PwENDXwN1fXVDC52FAscY0myIMkWRH58f7p8pcBSET6El2hInrT6+5macW7fifOddBnbuRA8Oomw2YlevJm7D9cRdd90lF20LBoIc3dbC7j/U4u7zUro8ndV3leBIvXD/kRCz3bTNBSemRlv/INVtZwPGGT591tDpGtG577BZKE6LY01JCsVpsRQnRzM3upNc1Yatv8YImO46OFwP3Q3nL0FsTzbCJLMCyu8YGTQJeWCJnpLPp7XmjOsM9X31Q6fNQqtoDh9pBpBuT6cwoZBbCm+hwFFgtGYcRWTHZWOewGm7i9YttDS18913cb7zLp5TpwCIys0l8SMfIW7D9cSsWIHJdumWndaa+spOtr9QRXeri6w5CWz+XAWZRTO2y1OIiJMAioCWXjc7qjvZWdPJzpouGrqGps6LtpgoSollbno8tyzIpCg11gibWC+JfcdRrfvhzGE4UQntx0euGGmxDbVa8q8d2YJJLACbY0o/l8vnGnm6bNhpM7d/aPaAswucVaRWcEfJHRQ6CsMtm9io2Eu8w+TwNjYy8P42Bra9z8COnQSdTrBYiFm+nPS77yZuw/VEFxWNev60M3V97HihiqaTPSSk27n1sQqKFqfK/GtCXIYE0DQ40zcYDpwdNZ3UdxqBk2CPYlVRMg9eW8C8zHiKUmPJdlgx9dZDa6Xxc/IwbK2EvmETjcdnGa2Y0g9BWtlQyMRljHnU2FgFggFaBlrOC5ravlraXG3hcmevmyl0FLI0Y+mIkMmIyZjWL+eAcwDX7l2h0NmGt74eAEt2Fo5bbyV23Tpi16654NLUl9Lb7mLnizVU7W3DHh/FdffOZcH6bMyypo0QoyIBNAXa+gbZEWrd7KzppLbDuHgy3mZhVVEKD15byOriZMpTojB1HIfW9+BkpRE0Z47A2etSlBlS50LBGiNwMhdCRoVxkeUU6/P2ndeSqe2tpaGvAW9waPRbfHQ8RY4iVmetDodMgaNgaLnmCNDBIINHjjKwbRsD77+P68AB8PtRdjuxK1eS9IlPELt2LdFFheMKQleflz2v1HFkaxMmi2L55kKu2ZRPtEydI8SYyG/MJKls7OU/9jSwvbqTmvZQ4FgtrCxK5oFV+awuTqHc3oO56nU4vQsOVkLHyaFrZKLjjYBZcl8obCogrXzSRpVdiNaa1oFWqnurqe6pDl8zU9dXN3I4s7KQG59LoaOQdTnrRrRmkm3JET/VpLXG19iIa9cuBrbvYGD7dgI9PQBY55eT8ulPE7t2Lfal12CKHn8/l3fQz6G3TrPv9Qb83iDz12ax4vYiYhOsk/RJhLi6SABNQDCoeet4Gz99r4ZdtV3ERJtZXZzClhV5rC5OYUFmHOaWfXDyl/DSq9AWupzKkQOZi6D8w0Nhk1gwJaPMwDht1uxsDgdNTW8NNT011PTW4PIP9T8l25IpdBSyIW+DETKhoMmNzyXKNLOGD/taWhjYtQvXrt0M7NqJv9mYa9aclkrc9dcTu24tsWvWYElJmfB7+b0BKt9tYt9r9Qw6fRQtTuXau0tIypz6/iohrmQSQOMw6Avwu32N/Oz9WmraB8hOsPHVzeXcuzIPBy6ofgs+eA1OvQ6uTuNUWsEa+NA3Ye4tkDJnSvpq/EE/p/tPU9NTMyJsantrR0yemR6TTnFCMXeX3m3M0pxYQnFCMUm2mTsljL+9nYFdu41Wzu5d+OobADAnJhKzciUxDz9M7OrVRBcXT1qLLOALcnRbM3v+WIer10ve/GRW3VFMRtHUDuYQ4mohATQGHU4P/76jnl/urKdrwEtFTgLf27KEzdkuoqr/CP/xKtRvN0am2ZOMQQJzb4aSG4zHk8QX8FHXV0d1bzW1PbXhsKnrqxtxgWZ2bDbFicWsylxlhExo8syxzgIQCf7ubly7P8C1aycDu3bjra4GwBQfT8yKFSTffz8xq1YZc6xNcssxGAhyfGcre/67jv6uQbLmJHDzIwvILp25AS3EbCQBNEqvVLbw//3nQVzeADeVp/PIdcWsspxCvf0YvBhacy+tHK79gtHKyV0xYtbl8fIFfJzqOcXRzqMc6TzCkY4jnOo5FQ4ahSI3PpeShBLW566nJLGEkoQSihKKiImKmfD7T4dwH87evbj37sO1fx/eKiNwVEwMMcuWkXjP3cSsXIVtfjnKPPnXAwH4fQGO72hl/+v19HUMkl7oYOMnysgtT4p4P5cQVyIJoMsIBjXfe/MU33vzFEvzE/nHjy5mju8EvP0oVL1hTPV/05Ow4C5jOPQE+II+qnuqOdJxJBw4J7tP4gvNWBAfHc/8lPl8cv4nmZc0j5LEEgodhREbbTZe2u9n8Nhx3Pv34dq7D9e+vQTaOwCjhWNfeg0Jt99BzMqV2CsWoqKmtv/J6/ZzeGsTB988javPS3qhg3UfK6VwkVzLI8RUkgC6BJfXz//8zUH+eLiVjyzN5ak1mug3/wxOvGLMKnDTk7DyzyB67J3R/qCf6p7qcNAc7TzKia4T4SHO8VFG2Hyi/BPMT53PguQF5MbnzsovxIBzAPfBA+HWjfvgIbTLGPwQlZ1N7OpriVm2FPs1S7GWzpn0U2oX4+73cvCt0xx+twmPy09uWRKbPjOfnHnS4hFiOkgAXUR7v4cHf76bE619fHNTBg90/jPq2RfBmgAb/xZWfXZMMwu4fC52t+5mR/MODnce5kTXifDAgNioWOanzOe+svtYkLqABSlG2IxmWYCZRgeDeOvqGTxcifvgIVz79+E5fgKCQTCZsJbNI/Huu43AWbqUqMzMaa9jT5uLQ283cuz9Zvz+IMVL0lh6cwEZhTK4QIjpJAF0AVprnvivSmranfzyYzms2fYw9DbC+r+Ga/981AMKTvefZmvjVt5rfI8PWj/AG/Rit9iZnzKfe+fdy/yU+SxIWUC+I392ho3W+M+cwV1ZyeChStyHKxk8fIRgv3EhrbLbsS9eTOpjj2FfuhT7ksVjnm1g0uoa1DQc66Ly7UbqD3diMinmrspg6c0FMpxaiAiRALqAPxxq4U9Hz/D0eitr3rkfPE745ItQcO0lt/MFfOxr22eETtN71PbWAlDoKOTesntZn7ueZenLiDLPrGtqRivQ04O78rDRuqk8jLvyULjvBosF29y5OG7bjL2iAtvCCqwlxShLZP+Led1+ju1oofKdRnrb3Ngd0ay4rZAF63PkAlIhIkwC6BxdA16+8fIR7s5s4+OH/5exeuen/9u4WPQCOtwdvNf4Hu81vcf25u0M+AaIMkWxInMFH5/7cdbnriffkT/Nn2Ligi4Xg8eODWvdHMbX0BB+Pbq4mLg1a7AtrMC+qAJrWRkm68z5Qu9uHaDynSaO72jB5wmQUeRg5WeKKFmajtky+1qbQlyJJIDO8b03ThI92M7/Nv8dyp4AD74IKSXnlWvoa+CZQ8/wh5o/ENAB0u3p3FJ4C+tz17M6a/WsGQINxkWeg8ePM3jsOJ7jxxg8dhxvXR2E1oqyZGVhr6gg8WMfNVo3CxZgjp951xL5PAGq97dxbFsLzad6MFkUpcsyqNiYK/07QsxAEkDDaA0vHWzm/6T+HnP/AHzyzfPCp76vnmcOPcN/1/w3FpOFLWVbuGvOXcxLmjfjR07pQABvfT2Dx47hCQXO4PHjBDo6wmWicnOxlZfhuP02bOXzsS+qwJKaGsFaX5rWmjO1fRzb3sKpPWfwDQZISLOz6s5i5q/NJsYxNWscCSEmTgJoGKfHR6n7BKuDfzQGG6SWhl+r7a3lmUPP8ErtK0Sborm//H4+veDTpMVM/czU4xEcGGDw5Ek8J06EguYYnhMn0YOhlVSjorCWziFu/XpsZWXYysuwzpuH2TE7WgquPi8ndrVybHsL3S0DWKJNzFmWTvmabLLmJMz4PwaEEBJAI/S4fDxp/X8QkwLXfxkAt9/NN3d+kz/U/AGr2conyz/JpxZ+ilT7zGgVnB2JNnj8OJ7jJxg8cRzPsePGmjehU2gmhwNbeTlJ996LtbwMW3k51qIi1ARmho4EvzdAXWUnJ3e3Ul/ZSTCoySw2ZiuYszydaJv8dxZiNpHf2GF8HjfLaILrvgW2BPq8fXzhzS9wsP0gD85/kE8t+BQp9onPrjxeQY8HT1XVUNAcP4HnxAkCvUNLb0fl5GAtL8Nx++3YQmFjycqatS2CQCBI47FuTn1whpqD7fgGA8Q4oll8Yx5la7JIzpIh1ELMVhJAw9iDTuNO+YfpdHfy2BuPUdVTxXfWf4cPFX5o2uqhtcbf3m6cPgu1bDwnT+CpqYVAAABls2GdO5f4D30Ia9k8bGVlWOfOnZGDA8ZKBzUt1T2c/KCN6n1tDDp9WGMszFmWTumKDHLmJmEyzc5AFWI0enp6eP755/n85z8/qvK1tbVs2bKFrq4uli5dyi9+8QuiL3CG47nnnuOb3/wmAH/7t3/LQw89NKbtJ5vSodM0AsqzY/UHX1lE/8Mv8OifHqV1oJXvbvwua3PWTtl76mAQb3U1g0ePMnj8BJ4Txxk8foJA17AF4bKysM2bh3XePGxl87DOKyO6IH/KJuWMBK017Q39nPrgDFV723B2e7BEmShanErpigzy56dgjpLh0+LqUFdXx+23387hw4dHVf7jH/8499xzD1u2bOGxxx5j8eLFfO5znxtRpquri+XLl7Nnzx6UUixbtoy9e/eSlJQ0qu0vRim1V2u9fMwfEgmgEZZnW/S7P/sqj+oqanpq+OGNP2RpxtJJfY9Aby/uQ4dw7z+A+8AB3IcOEXQaLS8VHY11zhysZWXhoLHNm4s5MXFS6zBTBANBmqt6qTnQTu2BdpzdHkwmRf6CZEpXZFC4KFX6dUTEPfn7Ixxt7pvUfc7PdvB3dyy46OtbtmzhpZdeYt68eWzatInvfOc7Fy2rtSYtLY3W1lYsFgs7duzgG9/4Bq+99tqIcr/61a945513+MlPfgLAZz/7WTZs2MCWLVtGtf3FTCSA5Ld7BM0Oh41DVYf4xrXfmHD4nG3duA6EwubAwfC6NphMWEtLcdx2G/bFi7FXLCS6qCjiMwdMNb83QMPRLmoPtlN3qJPBAR/mKBN55cmsvKOIokVp2OJm50wRQkyWb3/72xw+fJgDBw7Q39/PkiVLLlju+eefJz09ncTERCyh747c3FyamprOK9vU1EReXl748dlynZ2do9p+KlzZ33ZjpIF/bd1KoaOQO+fcOebtL9W6MSckYFuymITbb8O+ZAm2ikWY466ODvTBAR/1hzupOdBOw5FO/N4g0XYLhRUpFC9JI29+srR0xIx1qZbKdIiPj+fAgQMXfb29vf285y406OhCZ7uUUhd9fjrIb/0wvSYzXmcj/7zhn7GYRndofGfa6PvDH+j97z/gOXrMePJs62bzZuxLlmBfsoToosJZOxJtrLTWdDUPUH+kk4YjnbSc6iUY1MQkRFO2OoviJWlkz02UKXGEGIX+/n6uu+66C772/PPPU15eTk9PD36/H4vFQmNjI9nZ2eeVzc3N5Z133gk/bmxsZMOGDaSmpo5q+6kgATRMn9lMYUwGN+XfdMlyQZeL/jfeoPellxnYsQOCQWyLF5H2xb+86lo3Z3ncfhqPdYVCp4uBHmOpiZScWJZsyqNocRoZhQ6UjF4T4rLi4+PpD80qf7kWEMDGjRv57W9/y5YtW3juuee4887zz+DcfPPNPPHEE3R3dwPw+uuv89RTT6GUGtX2U0EGIQwTV2jXX/zl/+Af1v3Dea/pQADX7t30vvgSfX/6E9rlIionB8eH7yDhjg9jLS6KQI0jR2tNx2knDUc7qT/cSWtNHzqoibaZyStPJn9hCvnzU4hLmjkTlAoxm9x///0cOnSIW2+99ZKDEABqamrCw6ivueYafvnLX2K1WtmzZw8//vGPefbZZwH4+c9/zre+9S0AvvrVr/LpT3/6ktuPhoyCmyT2Irv+zRu/4Y6SO8LPBT0eOn/yE3p+9wL+M2cwxcXhuPUWEj78YezLlk3b6p0zQV+nm5ZTPTSe6KbhSBeuPmP11tS8OAoWpJC/IIWMYgdm89VzTIS42l01o+CUUsnAfwCFQB3wca119zll8oB/BzKBIPCM1vp7o32Pa9KvCd/3NjbS9JdfZPDoUeKuv56Ex79C3MaNmGy2CX+WmU5rTV+Hm6aTPTSf6qH5ZA/9XcY8ctZYC/nlyeQvSCFvfrKsqyOEGJdZ1QJSSv0j0KW1/rZS6nEgSWv9lXPKZAFZWut9Sql4YC9wl9b66OX2H1MUo/ur+zGbzDi3bqXpr78MwSDZT3+b+BtumJLPNFNorek546L5VE84dM7249jiosgpTSR7biLZpUmkZMdKX44QAriKWkDAncCG0P3ngHeAEQGktW4BWkL3+5VSx4Ac4LIBZEFhQtH+rz+g40c/wjpvHrnf/x7R+bNvQbnLCQSCdLcM0FLVa4TOqR7coVNqMY5osucmGqFTmkRSVsxVM4JPCDF9ZlsAZYQCBq11i1Iq/VKFlVKFwDXArtHs3IKJnt/8Jx0//CEJd99N5t99/Yo43RYManpaXbTV99FW309bfR8djU4CviAAcUlW8sqSyC5NJGduEgnpdgkcIcSUm3EBpJR6A6P/5lxfHeN+4oDfAX+ltb7oPBpKqUeBRwGS82Lp+OEPsS9bRta3/mFWfgnroKa33T0ibNpPO/F7jElMLVYz6fnxLLw+h/SCeDIKE3Ck2mblZxVCzG4zLoC01he9CEcpdUYplRVq/WQBbRcpF4URPv9Pa/3CZd7vGeAZgML0eO1vbyfnu/8yK76Qtdb0dw5ypq6P9vp+2hqMW+9gKGyiTKTmxTN/TRbpBfGkFThIzIiRmaSFEDPCjAugy3gZeAj4duj2pXMLKCM5fgYc01r/81h2HuP0E3f39cQsWzYZdZ1U7n4vnc0DdDU7jdsm4/7ZsDFZFKk5ccxdmUlaQTwZhQ6SMmMwyZBoIWadqVqOwWw2U1FRAUB+fj4vv/wyAD/4wQ/47ne/S3V1Ne3t7aSmTs+Cm7NtFFwK8BsgH2gAPqa17lJKZQPPaq03K6XWAe8BlRjDsAGe0Fq/crn9L7TZ9fb/egHHrbdO0Se4PK/bT1fLAJ1NTrqaB8Kh4+73hctYYy2kZMeRkh1LSm4c6QUOkrNjZWobIa4QU7EcA0BcXBzO0PyUw+3fv5+kpCQ2bNjAnj17xhRAciHqJFlos+v91VVE5eRM+Xv5fQG6W110DWvVdDY5cXZ5wmUsVjPJWbGkZMeSnB1LSnYcyTmxxDiiZ8UpQiGuCH98HForJ3efmRVw67cv+vJULMcAFw+gswoLC6c1gGbbKbgpFTQrLFMwCZ/X7ae9oZ+2hn7aQz+9bS7OZr/JrEjKjCWrJJGU9bEkh1o38ck2ud5GiKvQVCzHADA4OMjy5cuxWCw8/vjj3HXXXVP0CUZHAmiYoMU04ZaF1+2n/XQ/bfVDYdNzxhV+PS7JSlp+PHOWpYdbNQkZdpm+RoiZ6hItlekwWcsxADQ0NJCdnU1NTQ033HADFRUVlJSUTFZVx0wCaBg9xtbG2bBpbxgKnAuFzbxVGaTlO0jLjyfGMfXrrAshrhyTtRwDEH6+uLiYDRs2sH//fgmgmUKbLx1AWmvO1PZxbFszzVW9EjZCiCkxFcsxdHd3ExMTg9VqpaOjg23btvHlL395Kqo/ajIIYZjStAR9qr33vOcHB3yc2NXK0feb6WoewGI1k1eWRFp+POkFEjZCiMk32csxbN++nc9+9rOYTCaCwSB/9Vd/xcMPPwzA97//ff7xH/+R1tZW0tPT2bx5c3gJh8uRUXCTZE5Goq460wMYrZ2Wql6OvN9E9b52Ar4g6QXxzF+XTemKDFlCWgghkFFwk8g4BVdzoJ2dL1bT3eoi2mamfE0W89dlk5YXH+H6CSHElUMCaDilaK7q4bWfHiYxI4YbHixnzrJ0oqzmSNdMCCGuOBJAw2hMvPqTSuJTbNz9P5dii42KdJWEEOKKJRefDOMPJhAMaG77/CIJHyGEmGISQMNobWHDA2UkZcZGuipCCHHFkwAaQVOwMCXSlRBCiKuCBNAwSgVkwIEQIuJ6enr40Y9+NOryP/jBD5gzZw5KKTo6Oi5a7rnnnqO0tJTS0lKee+65MW8/2SSAhlHKH+kqCCHEmANo7dq1vPHGGxQUFFy0TFdXF08++SS7du1i9+7dPPnkk3R3d496+6kgo+CGUSp4+UJCiKvK07uf5njX8UndZ1lyGV9Z+ZWLvv74449TXV3NkiVLLrscA8A111xz2fd87bXX2LRpE8nJyQBs2rSJV199lfvuu29U208FCaARZFYIIUTkjWU5hvnz549qn01NTeTl5YUfX2rZhukiATSCBJAQYqRLtVSmw2gmIx2NC027FumFLSWAhBBiBrvccgyjbQHl5ubyzjvvhB83NjayYcOGSajh+EkADaOUtICEEJE31uUYRuPmm2/miSeeCA88eP3113nqqacmvN+JkFFwQggxw6SkpLB27VoWLlzIX//1X1+2/Pe//31yc3NpbGxk0aJFPPLIIwDs2bMnfD85OZmvfe1rrFixghUrVvD1r389PCDhYttPNVmOYZjirGxd09Ic6WoIIcSsMZHlGKQFJIQQIiIkgIQQQkSEBJAQQoiIkAASQggRERJAQgghIkICSAghRERIAAkhxAwzVcsxmM1mlixZwpIlS/jwhz885u0nmwSQEELMMFOxHAOA3W7nwIEDHDhwgJdffnnM2082mYpHCCEuofVb38JzbHKXY7CWl5H5xBMXfX0qlmOYyu3HSwJICCFmmKlYjgFgcHCQ5cuXY7FYePzxx7nrrrsmp8LjJAEkhBCXcKmWynSYrMlIARoaGsjOzqampoYbbriBiooKSkpKJmXf4yEBJIQQM9hkLccAkJ2dDUBxcTEbNmxg//79EkBCCCGGTMVyDN3d3cTExGC1Wuno6GDbtm18+ctfnvB+J0JGwQkhxAwzFcsxHDt2jOXLl7N48WI2btzI448/Hm49yXIMo6CUSgb+AygE6oCPa627L1LWDOwBmrTWt49m/7IcgxBCjM3VtBzD48CbWutS4M3Q44v5InBsWmolhBBizGZbAN0JPBe6/xxw14UKKaVygduAZ6enWkIIIcZqtgVQhta6BSB0m36Rct8FvgwEL7dDpdSjSqk9Sqk9gUBg0ioqhBDi0mbcKDil1BtA5gVe+uoot78daNNa71VKbbhcea31M8AzYPQBjb6mQgghJmLGBZDW+qaLvaaUOqOUytJatyilsoC2CxRbC3xYKbUZsAEOpdQvtdafmKIqCyGEGIfZdgruZeCh0P2HgJfOLaC1/hutda7WuhDYArwl4SOEEDPPbAugbwOblFKngE2hxyilspVSr0S0ZkIIMUnGOhv2Aw88wLx581i4cCGf+cxn8Pl8Fyz33HPPUVpaSmlpKc8991z4+draWlatWkVpaSn33nsvXq93wp9hNGZVAGmtO7XWN2qtS0O3XaHnm7XWmy9Q/p3RXgMkhBAzxXgC6Pjx41RWVuJ2u3n22fMHAHd1dfHkk0+ya9cudu/ezZNPPkl3t3EZ5Ve+8hW+9KUvcerUKZKSkvjZz342aZ/lUmZcH5AQQswk7/3mJB2nnZO6z9S8OK77+NyLvj7W5Rg2bx76+3vlypU0NjaeV+a1115j06ZNJCcnA7Bp0yZeffVVtmzZwltvvcXzzz8PwEMPPcQ3vvENPve5z43no42JBJAQQsww412Owefz8Ytf/ILvfe9755VtamoiLy8v/Dg3N5empiY6OztJTEzEYrGMeH46SAAJIcQlXKqlMh3GMhnp5z//edavX3/B2bMvNO2aUuqiz08HCSAhhJjBRrscw5NPPkl7ezs/+clPLlg2NzeXd955J/y4sbGRDRs2kJqaSk9PD36/H4vFQmNjY3jZhqkmASSEEDPMWJdjePbZZ3nttdd48803MZkuPLbs5ptv5oknnggPPHj99dd56qmnUEqxceNGfvvb37Jlyxaee+457rzzzkn9PBczq0bBCSHE1WCsyzE89thjnDlzhmuvvZYlS5bw93//98DI5RiSk5P52te+xooVK1ixYgVf//rXwwMSnn76af75n/+ZOXPm0NnZycMPPzx1H26YWbUcw1ST5RiEEGJsrqblGIQQQlwhJICEEEJEhASQEEKIiJAAEkIIERESQEIIISJCAkgIIURESAAJIYSICAkgIYQQESEBJIQQIiJkJoRhlFL9wIlI12OGSAU6Il2JGUCOwxA5FkPkWAyZp7WOH8+GMhnpSCfGO6XElUYptUeOhRyH4eRYDJFjMUQptWe828opOCGEEBEhASSEECIiJIBGeibSFZhB5FgY5DgMkWMxRI7FkHEfCxmEIIQQIiKkBSSEECIiJICEEEJExFUXQEqpW5RSJ5RSVUqpxy/wulJKfT/0+iGl1NJI1HM6jOJYPBA6BoeUUtuVUosjUc/pcLljMazcCqVUQCn10ems33QazbFQSm1QSh1QSh1RSr073XWcLqP4HUlQSv1eKXUwdCw+HYl6TjWl1M+VUm1KqcMXeX1835ta66vmBzAD1UAxEA0cBOafU2Yz8EdAAauBXZGudwSPxRogKXT/1qv5WAwr9xbwCvDRSNc7gv8vEoGjQH7ocXqk6x3BY/EE8HTofhrQBURHuu5TcCzWA0uBwxd5fVzfm1dbC2glUKW1rtFae4FfA3eeU+ZO4N+1YSeQqJTKmu6KToPLHgut9XatdXfo4U4gd5rrOF1G8/8C4C+A3wFt01m5aTaaY3E/8ILWugFAa32lHo/RHAsNxCulFBCHEUD+6a3m1NNab8X4bBczru/Nqy2AcoDTwx43hp4ba5krwVg/58MYf+FciS57LJRSOcDdwI+nsV6RMJr/F3OBJKXUO0qpvUqpB6etdtNrNMfiB0A50AxUAl/UWgenp3ozyri+N6+2qXjUBZ47dxz6aMpcCUb9OZVSGzECaN2U1ihyRnMsvgt8RWsdMP7YvWKN5lhYgGXAjYAd2KGU2qm1PjnVlZtmozkWNwMHgBuAEuBPSqn3tNZ9U1y3mWZc35tXWwA1AnnDHudi/OUy1jJXglF9TqXUIuBZ4Fatdec01W26jeZYLAd+HQqfVGCzUsqvtX5xWmo4fUb7O9KhtR4ABpRSW4HFwJUWQKM5Fp8Gvq2NjpAqpVQtUAbsnp4qzhjj+t682k7BfQCUKqWKlFLRwBbg5XPKvAw8GBrVsRro1Vq3THdFp8Flj4VSKh94AfjkFfjX7XCXPRZa6yKtdaHWuhD4LfD5KzB8YHS/Iy8B1ymlLEqpGGAVcGya6zkdRnMsGjBagiilMoB5QM201nJmGNf35lXVAtJa+5VSXwBewxjh8nOt9RGl1GOh13+MMcJpM1AFuDD+wrnijPJYfB1IAX4U+svfr6/AGYBHeSyuCqM5FlrrY0qpV4FDQBB4Vmt9weG5s9ko/1/8L+D/KqUqMU5DfUVrfcUt06CU+hWwAUhVSjUCfwdEwcS+N2UqHiGEEBFxtZ2CE0IIMUNIAAkhhIgICSAhhBARIQEkhBAiIiSAhBBCRIQEkBBCiIiQABJiiiil5iqlXldK9SqlfqGUske6TkLMJBJAQkyB0JXzfwI2AQHgExhT9wshQiSAhJga64B8jDnBMjHmxbo/ojUSYoa5qqbiEWIaFYZud2qtvUqpj0SyMkLMRBJAQkyNhNBtN0Boka5po5SqAbZrrT8xne8rxFjIKTghpkZ86LZ/rBsqpeYrpZ5TSp1WSnmVUh2hhd++o5SyjWL7WIwWWOVY31uI6SQtICGmRlzo1jmWjZRSdwD/CbQC/4axymQ2sB74LPDlUexmAcbMzBJAYkaTABJiapxtAY06gJRSDuA5jOBYr7V2n/N6th7d9PULQreHRvveQkSCnIITYmqMpwV0I5AE/Nu54QOgtb7QirV3KaX2KKXcSqlKpdQNwEKgW2vdGCqTpZRyKaWeP2fblaHnfzmGOgoxaSSAhJga4+kDigndzh9NYaXU54H/wjhd90WMhdNeAG5i2Om30MqUPwbuVUrNDW1bgLGK5R7gM2OooxCTRgJIiKkxnhbQ20Af8OdKqZNKqX9SSm1SSkWdW1AptRD4LvAvWuvbtdbPaK3/P+CHwCLO7/95GhgEngid6vtD6L3u1lp7x/LBhJgsEkBCTI0x9wGFTrGtAv4dSAf+B/A6UKeUuuuc4l8F3MA3znn+7dDtiP4frfUZ4P8AD2Asn5wFbNZad462fkJMNgkgIabGuEbBaa2Pa60fAlIxZlP4EZABPK+UyoTwND93AM9rrfvO2YU5dHuhEXD/hPE7vxK4S2tdNZa6CTHZJICEmBrjvg4IQGvt11pv01r/OfALwA4sCb1cAsQCey+w6XJAA4cv8NpXMX7nzYQukBUikiSAhJga42oBXYQndHs2zGIuVEgpZQEeBeq01v3nvPZXwJ8DX8Lo+/nGJNRLiAmRABJiasQBXq21bzSFlVLrLrRcg1JqDvAxoAHYFXq6IXR7wznF/yfGBKgjTr+FLm79J+AprfV3ge8BH1FKLR7dRxFiaqjRXdcmhBgtpZQVY8RZl9Y6ZZTbvA/Mw5gF4SDGH4cLgYcwlnO4VWu9fVj5PwK3AM9izLi9EWPph1Tgm1rrr4XKLQW2Ygw8uFdrrZVSSUA98JbW+q4Jf2AhxkkCSIhJppRKBdqBeq114Si3uRv4CMYAgSwgCqOl8yrwHa316XPKp2Nc23MjRkC9i9FX9DuMoPmNUioXo9V0GtigtR4ctv23gL8BlmutL9SXJMSUkwASYpIppQqBWuCI1nphhKsjxIwlfUBCTAKlVIFS6kOh2aqTQ0+fO0RaCDGMBJAQk+MTGFPh3ACsDj1XE7nqCDHzySk4ISaBUqoMY/YBL0b/TTRwm9b6lYhWTIgZTFpAQkwCrfVx4GGMa3WagAclfIS4NGkBCSGEiAhpAQkhhIgICSAhhBARIQEkhBAiIiSAhBBCRIQl0hUQYrZTSl0xI3m01irSdRBXDxkFJ4QQIiLkFJwQM4RS6iallFZKPRLpuggxHSSAhJg5FoVuD0W0FkJMEwkgIWaOxUCQCy+nLcQVRwJIiJljMVCltXZFuiJCTAcJICFmAKVUFFCOsRqqEFcFCSAhZoYyjBm0pf9HXDUkgISYGRaHbqUFJK4aEkBCzAxnA0haQOKqIQEkxMywGOjRWtdHuiJCTBcJICFmhkVA5aUKKKX+VSn13LDHv1JK/XzKaybEFJEAEiLClFIZQAaX7//5FnCPUqpYKfVkaJvPTnX9hJgqMhmpEJF3tv+nUCn1+AVe/63Wukpr3aKU+inwAmAF1mitfdNWSyEmmQSQEJF3NoBuD/2c6/fD7h8AvgR8SGvdPcX1EmJKyWzYQswSSqnlwEvAViCgtf5EhKskxIRIAAkxCyil8oH3gc8A+4EaYJXW+nhEKybEBEgACTHDKaUcwDbg+1rrn4ae+w6Qq7W+L6KVE2ICJICEEEJEhAzDFkIIERESQEIIISJCAkgIIURESAAJIYSICAkgIYQQESEBJIQQIiIkgIQQQkSEBJAQQoiIkAASQggREf8/nzyrUPgDIXkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(0,21,5):\n",
    "  plt.plot(S_ave[i,:],y,label='t=%4.2f' %t[i])\n",
    "\n",
    "plt.ylim([-0.5,0.5])\n",
    "plt.xlim([0,1])\n",
    "plt.xlabel(r'$\\frac{\\int \\ S dx}{L_x}$',fontsize=24)\n",
    "plt.ylabel(r'$y$',fontsize=24)\n",
    "plt.legend(loc='lower right').draw_frame(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:dedalus]",
   "language": "python",
   "name": "conda-env-dedalus-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
